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TECHNICAL PUBLICATION 


PROBABILITY AND STATISTICS IN AEROSPACE ENGINEERING 

I. INTRODUCTION 


A. Preliminary Remarks 

Statistics is the science of the collection, organization, analysis, and interpretation of numerical 
data, especially the analysis of population characteristics by inference from sampling. In engineering work 
this includes such different tasks as predicting the reliability of space launch vehicles and subsystems, life- 
time analysis of spacecraft system components, failure analysis, and tolerance limits. 

A common engineering definition of statistics states that statistics is the science of guiding 
decisions in the face of uncertainties. An earlier definition was statistics is the science of making decisions 
in the face of uncertainties, but the verb making has been moderated to guiding. 

Statistical procedures can vary from the drawing and assessment of a few simple graphs to carry- 
ing out very complex mathematical analysis with the use of computers; in any application, however, there 
is the essential underlying influence of “chance.” Whether some natural phenomenon is being observed or 
a scientific experiment is being carried out, the analysis will be statistical if it is impossible to predict the 
data exactly with certainty. 

The theory of probability had, strangely enough, a clearly recognizable and rather definitive start. It 
occurred in France in 1654. The French nobleman Chevalier de Mere had reasoned falsely that the prob- 
ability of getting at least one six with 4 throws of a single die was the same as the probability of getting at 
least one “double six” in 24 throws of a pair of dice. This misconception gave rise to a correspondence 
between the French mathematician Blaise Pascal (1623-1662) and his mathematician friend Pierre Fermat 
(1601-1665) to whom he wrote: “Monsieur le Chevalier de Mere is very bright, but he is not a mathe- 
matician, and that, as you know, is a very serious defect.” 

B. Statistical Potpourri 

This section is a collection of aphorisms concerning the nature and concepts of probability and 
statistics. Some are serious, while others are on the lighter side. 

“The theory of probability is at bottom only common sense reduced to calculation; it makes us 
appreciate with exactitude what reasonable minds feel by a sort of instinct, often without being able to 
account for it. It is remarkable that this science, which originated in the consideration of games of chance, 
should have become the most important object of human knowledge.” (P.S. Laplace, 1749-1827) 

“Statistical thinking will one day be as necessary for efficient citizenship as the ability to read 
and write.” (H.G. Wells, 1946) 



From a file in the NASA archives on “Humor and Satire:” Statistics is a highly logical and 
precise method of saying a half-truth inaccurately . 

A statistician is a person who constitutionally cannot make up his mind about anything and under 
pressure can produce any conclusion you desire from the data. 

There are three kinds of lies: white lies, which are justifiable; common lies, which have no 
justification; and statistics. 

From a NASA handbook on shuttle launch loads: “The total load will be obtained in a rational 
manner or by statistical analysis.” 

Lotteries hold a great fascination for statisticians, because they cannot figure out why people play 
them, given the odds. 

There is no such thing as a good statistical analysis of data about which we know absolutely 
nothing. 

Real-world statistical problems are almost never as clear-cut and well packaged as they appear 
in textbooks. 

Statistics is no substitute for good judgment. 

The probability of an event depends on our state of knowledge (information) and not on the state of 
the real world. Corollary: There is no such thing as the “intrinsic” probability of an event. 

C. Measurement Scales 

The types of measurements are usually called measurement scales. There exist four kinds 
of scales. The Ust proceeds from the “weakest” to the “strongest” scale and gives an example of each: 

• Nominal Scale: Red, Green, Blue 

• Ordinal Scale: First, Second, Third 

• Interval Scale: Temperature 

• Ratio Scale: Length. 

Most of the nonparametric (distribution-free) statistical methods work with interval or ratio scales. 
In fact, all statistical methods requiring only a weaker scale may also be used with a stronger scale. 

D. Probability and Set Theory 

The formulation of modem probability theory is based upon a few fundamental concepts of set 
theory. However, in probability theory these concepts are expressed in a language particularly adapted to 
probability terminology. In order to relate the notation commonly used in probability theory to that of set 
theory, we first present a juxtaposition of corresponding terms, shown in table 1. 
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Table 1.— Set theory versus probability theory terms. 


Set Vocabulary 

Probability Vocabulary 

(1) Element 

Outcome (E) 

(Sample Point, Elementary Event) 

(2) Subset 

Event (A) 

(3) Universal Set 

Sample Space ( S) 

(4) Empty Set 

Null Event (<P) 

(5) Disjoint 

Mutually Exclusive 

(6) Union AkjB 

“OR” Probability 

(7) Intersection AnB 

“AND” Probability 


The probability theory terms in table 1 are defined as follows: 

(1) An outcome E is defined as each possible result of an actual or conceptual experiment. Each 
experiment terminates with an outcome. An outcome is sometimes called a sample point or elementary 
event. 


(2) An event A is defined as a set of outcomes. One declares that an event A has occurred if an 
outcome E of an experiment belongs to an element of A. 

(3) The sample space S is defined as the set of all possible outcomes. It is also called the certain 

event. 

(4) The null event 0 is defined as the set consisting of no outcomes. It is also called the 
impossible event. 

(5) Two events A and B are called mutually exclusive if they have no common element. Note that 
outcomes are by definition mutually exclusive. 

(6) The union of events A and B is the event that occurs if A occurs or/and B occurs. 

(7) The intersection of events A and B is the event that A occurs and B occurs. 

Two more definitions are used in probability theory with notations that are identical to that of set 

theory: 

(8) The complement of an event A, written as A, A c ,or A', is the event that occurs if A does not 

occur. 


(9) The difference of events A and B, written as A-B, is the event that occurs if A occurs but B 
does not occur: (A-B)=An B ' . 

EXAMPLE: Toss a die and observe the number that appears facing up. Let A be the event that an even 
number occurs, and B the event that a prime number occurs. Then we have in table 2: 
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Table 2 . — Examples of set theory. 


Sample Space 

S={1 ,2,3,4,5,61 

Outcome 

£=m,{2),(3}.{4},{5},{6} 

Event A 

A={2,4,6} 

Event B 

£={2,3,5} 

Union 

AuB={ 2,3, 4,5,6} 

Intersection 

Anfl={ 2} 

Complement 

^'={1 ,3,5} 

Difference 

A- £={4,6} 


1. Venn Diagrams 

When considering operations on events it is often helpful to represent their relationships by 
so-called Venn Diagrams, named after the English logician John Venn (1834-1923). 

2. Principle of Duality (De Morgan’s Law) 

The Principle of Duality is also known as De Morgan’s Law after the English mathematician 
(1871). Any result involving sets is also true if we replace unions by intersections, intersections by 
unions, and sets by their complements. For example, (A'uB)'=A'nB'. 
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II. PROBABILITY 


A. Definitions of Probability 


1. Classical (a Priori) Definition 

The classical (a priori) definition of probability theory was introduced by the French mathe- 
matician Pierre Simon Laplace in 1812. He defined the probability of an event A as the ratio of the 
favorable outcomes to the total number of possible outcomes, provided they are equally likely (probable): 

P(A)=jj (1) 

where n=number of favorable outcomes and jV=number of possible outcomes. 


2. Empirical (a Posteriori) Definition 

The empirical (a posteriori) definition was introduced by the German mathematician Richard V. 
Mises in 1936. In this definition, an experiment is repeated M times and if the event A occurs m(A) times, 
then the probability of the event is defined as: 


P( A)— lim 

M—>oo 


m(A) 

M 


( 2 ) 


Empirical Frequency. This definition of probability is sometimes referred to as the relative 
frequency. Both the classical and the empirical definitions have serious difficulties. The classical definition 
is clearly circular because we are essentially defining probability in terms of itself. The empirical definition 
is unsatisfactory because it requires the number of observations to be infinite; however, it is quite useful in 
practice and has considerable intuitive appeal. Because of these difficulties, statisticians now prefer the 
axiomatic approach based on set theory. 


3. Axiomatic Definition 

The axiomatic definition was introduced by the Russian mathematician A.N. Kolmogorov in 

1933: 


• Axiom 1: P(A )> 0 

• Axiom 2: P(S)= 1 

• Axiom 3: P(AuB)=P(A)+P(£) if AnB=0 . 

It follows from these axioms that for any event A, then: 

0<P(A)<1 . (3) 
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Probabilities and Odds. If the probability of event A is p, then the odds that it will occur are 
given by the ratio of p to 1 -p. Odds are usually given as a ratio of two positive integers having no com- 
mon factors. If an event is more likely to not occur than to occur, it is customary to give the odds that it 
will not occur rather than the odds that it will occur. 

EXAMPLE: 

• Probability: P=A/B where A and B are any two positive numbers and A<B . 

Odds: A: ( B-A ). 

If the probability of an event is p= 3/4, we say that the odds are 3:1 in its favor. 

• Given the odds are A to B, then the probability is A/(A+B). 

Criticality Number. For high reliability systems it is often preferable to work with the probability 
of failure multiplied by 10 6 . This is called the criticality number. For instance, if the system has a 
probability of success of P=0.9999, then the criticality number is C=100. 


B. Combinatorial Analysis (Counting Techniques) 

In obtaining probabilities using the classical definition, the enumeration of outcomes often becomes 
practically impossible. In such cases use is made of combinatorial analysis, which is a sophisticated way 
of counting. 


1. Permutations 

A permutation is an ordered selection of k objects from a set S having n elements. 
• Permutations without repetition: 


P 0 (n, k)= n P k =nx(n-l)x(n-2)...(n-k+l)=j^jj 


( 4 ) 


• Permutations with repetition: 


Pl(n,k)=n^ . 


( 5 ) 


. Combinations 

A combination is an unordered selection of k objects from a set S having n elements. 

• Combinations without repetition: 

C Q (n,k)= n C^=^j= k f (n-k) f ( ca ^ ec * “binomial coefficient”) . (6) 
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Combinations with repetition: 


c 1 K*j=C ,+ *- 1 ) 


(ti+k-1)! 

k!(n-\)l ' 


EXAMPLES: Selection of two letters from {a, b, c}: 


(7) 


(1) Po(n,k)=Po (3, 2)=3x2=6 Without repetition 

ab, ac, ba, be, ca, cb 

(2) P\(n, k)=Pi(3, 2)=3x3=9 With repetition 

aa, ab, ac, ba, bb, be, ca, cb, cc 

(3) C 0 (n, k)=Co (3, 2)= =3 Without repetition 

ab, ac, be 

(4) Ci(«,ifc)=Ci (3,2)=-gf=6 With repetition 

aa, bb, cc, ab, ac, be 

PROBLEM: Baskin-Robbins ice-cream parlors advertise 31 different flavors of ice cream. What is the 
number of possible triple-scoop cones without repetition, depending on whether we are interested in how 
the flavors are arranged or not? 

SOLUTION: 3 iP 3 =26,970 and 3 iC 3 =4,495. 


3. Permutations of a Partitioned Set 


Suppose a set consists of n elements of which n\ are of one kind, «2 are of a second kind, . . .n^ 
are of a k^ type. Here, of course, n=n\+n 2 +...ni i . Then the number of permutations is: 


P 2 (n,n k )= 


n[ 

n l !n 2 !n 3 !...n k ! 


( 8 ) 


An excellent reference for combinatorial methods is M. Hall’s book Combinatorial Analysis. 

PROBLEM: Poker is a game played with a deck of 52 cards consisting of four suits (spades, clubs, 
hearts, and diamonds) each of which contains 13 cards (denominations 2, 3, 4, 5, 6, 7, 8, 9, 10, J, Q, K, 
and A.) When considered sequentially, the A may be taken to be 1 or A but not both; that is, 10, J, Q, K, 
A is a five-card sequence called a “straight,” as is A, 2, 3, 4, 5; but Q, K, A, 2, 3 is not sequential, that is, 
not a “straight.” 

A poker hand consists of five cards chosen at random. A winning poker hand is the one with 
a higher “rank” than all the other hands. 
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A “flush” is a five-card hand all of the same suit. A “pair” consists of two, and only two, cards of 
the same kind, for example ( J s , J c ). “Three -of- a-kind” and “four-of-a-kind” are defined similarly. A “full 
house” is a five-card hand consisting of a “pair” and “three-of-a-kind.” The ranks of the various hands are 
as follows with the highest rank first: 

• Royal flush (10, J, Q, K, A of one suit) 

• Straight flush (consecutive sequence of one suit that is not a royal flush) 

• Four-of-a-kind 

• Full house 

• Flush (not a straight flush) 

• Straight 

• Three-of-a-kind (not a full house) 

• Two pairs (not four of a kind) 

• One pair 

• No pair (“bust”). 

(1) Show that the number of possible poker hands is 2,598,960. 

(2) Show that the number of possible ways to deal the various hands are: 

(a) 4 (b) 36 (c) 624 (d) 3,744 

(f) 10,200 (g) 54,912 (h) 123,552 (i) 1,098,240 

SOLUTIONS: Poker is a card game with a deck of 52 cards consisting of four suits (spades, clubs, 
hearts, and diamonds). Each suit contains 13 denominations (2, 3, 4, 5, 6, 7, 8, 9, 10, J, Q, K, and A). 
A poker hand has five cards and the players bet on the ranks of their hands. The number of possible ways 
to obtain each of these ranks can be determined by combinatorial analysis as follows: 

(a) Royal Flush. This is the hand consisting of the 10, J, Q, K, and A of one suit. There are four 
of these, one for each suit, and hence, N= 4. 

(b) Straight Flush. All five cards are of the same suit and in sequence, such as the 6, 7, 8, 9, and 
10 of diamonds. Their total number is 10 for each suit. However, we have to exclude the four royal 

flushes contained in this set. Therefore, the total number of straight flushes is M=10x4-4=36. 

(c) Four-of-a-Kind. This hand contains four cards of the same denomination such as four aces or 
four sixes and then a fifth unmatched card. The total number of possible ways to choose this rank is 
obtained by: 

(1) Choosing the denomination, 13 ways. 

(2) Choosing the suit, 4 ways. 

(3) Choosing the remaining unmatched card, 12 ways. 

The result is: A=13x4x 12=624. 

(d) Full House. This hand consists of three cards of one denomination and two cards of another, 
as 8-8-8-K-K. The total number of possible ways is given by the following sequence of selections: 

(1) Choosing denomination of the first triplet of cards, 13 ways. 


(e) 5,108 
(j) 1,302,540 
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(2) Selecting three out of the four suits for this triplet, ( 3 ) =4 ways. 

(3) Choosing the denomination of the second doublet of cards, 12 ways. 

(4) Selecting two out of the four suits for this doublet, wa y s - 

The result is then N=13x4xl2x6=3,744. 


(e) Flush. This hand contains five cards of the same suit, but not all in sequence. To obtain the 
number of possible ways, select 5 out of 13 denominations, | 1 5 3 j=l,287 ways, and select one of the four 
suits, 4 ways for a total of N=l, 287x4 =5,148. Here we have to consider again that this number contains 
the number of straight flushes and royal flushes which have to be subtracted from it. The result is 

A=5, 148-36^1=5, 108. 

(f) Straight. This hand contains a five-card sequence as defined above. We observe that there are 
10 possible five-card sequences. Each face value in this sequence can come from any of the four denomi- 
nations, which results in 4 5 different ways of creating a particular five-card sequence. The result is: 
N=10x4 5 = 10,240. Again, it must be noted that this number contains the number of straight flushes and 
royal flushes which have to be subtracted for the final answer: N= 10,240- 36 -4= 10,200. 

(g) Three-of-a-Kind. This hand contains three cards of one denomination and two different cards, 
each of different denominations. The total number of ways is obtained by: 

(1) Choose the denomination of the three cards, 13 ways. 

(2) Select three of the four suits in 4 ways. 

(3) Select 2 of the remaining 12 denominations |^j =66 ways. 

(4) Each of the two remaining cards can have any of the four denominations for 4x4=16. 

The total number for this rank is, therefore, N= 1 3 x4x66x 1 6=54,9 1 2. 

(h) Two Pairs. To obtain the number of possible ways for this rank, we take the following steps: 

(1) Select the denomination of the two pairs in ^ )=78 ways. 

(2) Select the two suits for each pair in ^j=36 ways. 

(3) Select the denomination of the remaining card. There are 1 1 face values left. 

(4) The remaining card can have any of the four suits. 
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The total number is, therefore, 7^=78x36x1 1x4=123,552. 

(i) One Pair. The number of possible ways for this rank is obtained according to the following steps: 

( 1 ) Select denomination of the pair in 1 3 ways. 

(2) Select suit in wa y s - 

(3) Select denomination of the other three cards from the remaining 12 denominations 
in | j =220 ways. 

(4) Each of these three cards can have any suit, resulting in 4 3 =64 ways. 

The total number is then 7V= 13x6x220x64=1 ,098,240. 

(j) No Pair. The number of ways for this “bust” rank is obtained according to the following steps: 

(1) Select five cards from 13 denominations as ^j=l,287. 

(2) Each card can have any suit, giving 4 5 = 1,024. 

The result is 7V=1, 287x1, 024=1, 317, 888. Again, we note that this number contains the number 
of royal flushes, straight flushes, flushes, and straights. Thus, we obtain as the answer: 

7V=1,3 17,888-4-36-5,108-10,200=1,302,540 . 

QUESTION: The Florida State Lottery requires the player to choose 6 numbers without replacement from 
a possible 49 numbers ranging from 1 to 49. What is the probability of choosing a winning set of num- 
bers? Why do people think a lottery with the numbers 2, 13, 17, 20, 29, 36 is better than the one with the 
numbers 1, 2, 3, 4, 5, 6? (Hint: use hypergeometric distribution.) 

C. Basic Laws of Probability 

1. Addition Law (“OR” Law; “AND/OR”) 

The Addition Law of probability is expressed as: 

P(AuB)=P(A)+P(B)-P(AnB) . (9) 

The Venn diagram in figure 1 is helpful in understanding this probability law. A formal proof can 
be found in any standard textbook. In Venn diagrams, the universal set S is depicted by a rectangle and the 
sets under consideration by closed contours such as circles and ellipses. 
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GENERAL RULE: 



n= 3: P(AuB^C)=P(A)+P(B)+P(C)-P(Ar^B)-P(Ar^C)-P(BnC)+P(AnBnC) (10) 

n arbitrary (obtainable by mathematical induction): 

k k k 

P(A l uA 2 v...vA k )=? l P(A i h I P(A i nA j )+ £ P(A j nA j nA r ) 

1 i<j=2 i<j<r = 3 ' 

Af-l)** 1 P(A l nA 2 nA 3 n...nA k ) . (11) 


2. Conditional Probability 

Since the choice of sample space is not always self-evident, it is often necessary to use the symbol 
P(A I B) to denote the conditional probability of event A relative to the sample space B, or the probability 
of A given B. Assuming equal probability for the outcomes in A and B, we can derive the relationship 
shown in figure 2. 



Figure 2. — Conditional probability. 


Given the number of outcomes in sample space B as N(B), the number of outcomes in sample 
space Ar\B as N(Ar\B), and the number of outcomes in the sample space S as N(S), we obtain P(A I B) 
using the classical definition of probability: 


P,A I B> _ A^(AnR) _ N(Angy^) 

N(B ) N(B)/N(S) 


( 12 ) 


The second term, on the right-hand side, was obtained by dividing the numerator and denominator 
by N(S). The sample space B is called the reduced sample space. 





We can now write the above equation in terms of two probabilities defined on the total sample 
space S: 


P(A I B)= 


P(AnB) 

P(B) 


(13) 


Generalizing from this example, we introduce the following formal definition of conditional 
probability: If A and B are any two events in a sample space S and P(B)*0, the conditional probability 
of A given B is: 


P(A\B) = 


P(AnB) 
Pi B) 


(14) 


3. Multiplication Rule (“AND” Law) 

Multiplying the above expression of the conditional probability by P(B), we obtain the following 
multiplication rule: 


P(A nB) = P{B) P(A I B) . (15) 

The second rule is obtained by interchanging letters A and B. This rule can be easily generalized 
to more than two events; for instance, for three events A, B, and C we have: 

P(AnBnC) = P(A)P(B\A)P(C\AnB) . (16) 

If A and B are two events, we say that A is independent of B if and only if P(A\B)=P(A). In other 
words, the occurrence of event B does not affect the probability of event A. It can be shown that B is inde- 
pendent of A whenever A is independent of B. Therefore, we can simply state that A and B are 
independent events. 

EXAMPLE: Two people often decide who should pay for the coffee by flipping a coin. To eliminate the 
problem of a biased coin, the mathematician John V. Neumann (1903-1957) devised a way to make the 
odds more “even” using the multiplication law. 

The coin is flipped twice. If it comes up heads both times or tails both times, the process is 
repeated. If it comes up heads-tail, the first person wins, and if it comes up tail-heads the second person 
wins. The probabilities of these outcomes are the same even if the coin is biased. 

For example, if the probability of heads is P(H)= 0.6 and that of tails is P(7)=0.4, the intersection 
of the two events P(//n7')=P(7’n//)=0.6x0.4=0.24. 

4. Event-Composition Method 

This approach calculates the probability of an event A by expressing it as a composition (unions 
and/or intersections) of two or more other events. 

PROBLEM 1: Use the addition and multiplication laws of probability to simplify the expression: 

P[(BnC)v(DnE)] . (17) 
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SOLUTION: Applying the addition law, we obtain: 


P[(BnC)u(DnE)]=P(BnC)+P(DnE)-P(BnCnDnE) . (18) 

Observe that the events on the right are intersections and this calls for the application of the multiplication 
law. Therefore, 


P[(BnC)u(DnE)]=P(B) P(C\B)+P(D) P(E\D)-P(B) P{C\B) P(D\BnC) P(E\BnCnD ) . (19) 

It is frequently desirable to form compositions of mutually exclusive or independent events, 
because they simplify the addition and multiplication laws. 

PROBLEM 2: It is known that a patient will respond to a treatment of a particular disease with a probabil- 
ity of 0.9. If three patients are treated in an independent manner, determine the probability that at least one 
will respond. 

SOLUTION: Define the events: 

A=At least one patient will respond. 

A ( =/th patient will respond (i=l, 2, 3). 


The event A=AiuA 2 uA 3 . 


Now we observe by the law of duality that the complementary event A' is Aj'nA^nAj and that 
S-AuA' . Then, because P(S)= 1 and the independence of the events A,- W e have: 


P(A)^l-P(A') = l-P(A{)xP(A^)xP(A^) 
P(A)=l-0.1x0.1x0.1=0.999 . 


( 20 ) 


This result is of importance because it is often easier to find the probability of the complementary event A' 
than of the event A itself. This is always the case for problems of the “at-least-one” type, as is this one. 

PROBLEM 3: Birthday Problems: (a) What is the probability that in a randomly selected group of 
n persons, two or more of them will share a birthday? 

SOLUTION: In solving this problem, we assume that all birthdays are equally probable (uniform distri- 
bution). We also discount leap years. These assumptions are, of course, not quite realistic. Again, it is 
advantageous to first find the complementary event that all persons have different birthdays. 

The first of the n persons has, of course, some birthday with probability 365/365=1. Then, if the 
second person is to have a different birthday, it must occur on one of the other 364 days. Thus the proba- 
bility that the second person has a birthday different from the first is 364/365. Similarly the probability that 
the third person has a different birthday from the first two is 363/365 and so on. 

The probability of the complementary event A' is, therefore: 

P( A')= (365/365)x(364/365 )x. . .( (365-n+l )/365 ) . (21) 
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The desired probability of the event A is, then: 

P(A)=\-P(A') . (22) 

For «=23: P(A)=0.5073 and for n= 40: P(A)=0.891. 

(b) What is the probability that in a randomly selected group of n persons, at least one of them will 
share a birthday with you ? 

SOLUTION: The probability that the second person has a birthday different from you is, of course, the 
same as above, namely 364/365. However, the probability that the third person has a different birthday 
from yours is, again, 364/365 and so on. 

The probability of the complementary event A' is, therefore: 

P(A') = (364/365 )” -1 . (23) 

The desired probability of event A is, then, again: 

P(A)=l-P(A') . (24) 

For n=23\ P(A)= 0.058 and for n= 40: P(A)=0.101. 

PROBLEM 4: Three cards are drawn from a deck of 52 cards. Find the probability that two are jacks and 
one is a king. 

SOLUTION: p=(3 !/2 ! )x(4/52x3/5 1 )x(4/50)=6/5525 . 

5. Total Probability Rule 

Let B \ , B 2 , ■■■, Bk form a partition of the sample space S (see fig. 3). That is, we have 


BiC\Bj=0 for all tej 


and 


B\uB2U...uB/ ( =S . 

The events 5, are mutually exclusive and exhaustive (see fig. 3). 


(25) 

(26) 
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Figure 3. — Partitioned sample space. 

The total probability for any event A in S is: 

k k 

P(A)='ZP(AnB i )=J t P(B i )xP(A\B i ) . (27) 

(=1 i=l 

6. Bayes’ Rule 

Bayes’ Rule was published in 1763 by Thomas Bayes. Bayes’ formula (see fig. 4) finds the 
probability that the event A was “caused” by the events (inverse probabilities). 



Figure 4. — Bayes’ Rule. 


Let Bj form a partition of S: 


P(B.\A)= 


P(AnB t ) _ P(A I BjjxPfBi) 
P(A) ~ P(A) 


Substituting in the denominator: 


( 28 ) 
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P(A)=ZP(AnB i )=ZP(A\B i )xP(B i ) , 


(29) 


we obtain Bayes’ formula as: 


P(A\B ; )P(B ; ) 
P(B i \A)=-£- l 

J d P(A\B i )P(B i ) 


(30) 


The unconditional probabilities P(B,) are called the “prior” probabilities. The conditional 
probabilities P(B, A) are called the “posterior” probabilities. 

PROBLEM 1 : Suppose a person with AIDS is given an HIV test and that the probability of a positive 
result is 0.95. Furthermore, if a person without AIDS is given an HIV test, the probability of an incorrect 
diagnosis is 0.01. If 0.5 percent of the population has AIDS, what is the probability of a healthy person 
being falsely diagnosed as being HTV positive? 

SOLUTION: Define the events: 


A=Person has AIDS 
B=HIV test is positive 


P(B I A)= 0.95 P(B I A') = 0.01 

P(A)=0.005 P(A') = 0.995 


P(B\A')xP(A') (0.01)(0.995) 

' ' P( B I A' )x P( A' )+ P( B \ A)xP( A) (0.01)(0.995)+(0.95)(0.005) 


= 0.6768 . 


(31) 


MESSAGE: False positive tests are more probable than the true positive tests when the overall population 
has a low incidence of the disease. This is called the false-positive paradox. 


NOTE 1: In the medical profession P(A) is called the base rate and the event B I A' is called a false positive 
test (later to be defined as a type I error) and B' I A is called a false negative test (type II error). 

NOTE 2: Often it is said that a test proves to be 95 percent reliable or accurate. This statement means that 
P(B\A)=P(B' I A')=0.95, but it is more precise in this case to call the test 95 percent accurate in both 
directions. 


AUXILIARY QUESTION: What is the probability of a person with AIDS being incorrectly diagnosed as 
not being HTV positive (type II error)? 


P(A I B')= 


P(B'\A)xP(A) 


(0.05)(0.005) 


P(B' I A)xP(A)+P(B' I A')xP(A') (0.05)(0.005)+(0.99)(U995) 


= 2.537x 1CT 4 . (32) 


SOLUTION: False negative tests are highly improbable. 

PROBLEM 2: Suppose you are on a game show (Let’s Make a Deal), and you are given the choice of 
three doors. Behind one door is a car; behind the other two, goats. You pick a door, say #1, and the host 
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(Monty Hall), who knows what is behind the doors, opens another door, say #3, which has a goat. He 
then says to you, “Do you want to pick door #2?” Is it to your advantage to switch your choice? 

SOLUTION: Define events: 

• Di = Car is behind door D/(*=l, 2, 3) 

• // 3 =Host opens door #3. 

Define conditional probabilities: 

• P(Z> 2 1 // 3 )=Probability that car is behind door #2, given host has opened door #3 

• P (/?3 I D,)=Probability that host opens door #3, given that car is behind door #1 . 

Apply Bayes’ formula: 


P(D 2 \H 3 ) = 


P(H 3 \D 2 )xP(D 2 ) 

P(W 3 I £>j jxPfDj J+P(7r 3 I £> 2 )XP(Z> 2 J+P(jr 3 I D 3 )XP(D 3 ) • 


Prior probabilities that car is behind door #1 : 


P(D l )=P(D 2 )=P(D 3 )=V3 . 

Conditional probabilities: 

P(Hi \D\)=l/2 (some set this to unknown q), P(H 3 ID2)=1, P{Ht, l£> 3)=0 . 
Posterior probability: 


(33) 


(34) 

(35) 


P(D 2 \H 3 ) = 




(l/2)x(l/3)x(l)x(l/3)+(0)x(l/3) 


=2/3 


(36) 


Therefore, it is of advantage to switch. 

PROBLEM 3: We consider 10 successive coin tosses. If the coin is “fair” and the events are assumed to be 
independent, then the probability of obtaining 10 heads in a row is clearly given as Pj o=( 1/2) 1 °= 1/1 024. 
What is the probability that the 1 1th toss will be a head? 

SOLUTION: With the above assumptions, the answer is obviously one-half. Sometimes it is said that the 
coin has no “memory.” (There are some gamblers who will bet on tails because of the law of averages, 
thinking it acts like a rubber band: the greater the deviation, the greater the restoring force towards the 
mean. This is also known as the “gambler’s fallacy”.) However, it is more natural to think that the longer 
the run of heads lasts, the more likely it is that our assumption of the coin being “fair” is not true. Let us 
expand this train of thought by defining the following situation: 


Definition of events: 

F=Coin is unbiased 
B=Coin is biased 
H= Next toss will be heads. 
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Definition of probabilities: 


P(B)= 1 -P(F)=0. 1 0 
P(H I B)=0.70 . 


F(F)= 0.90 
P(H I F)=0.50 

Applying the total probability rule we obtain: 


(37) 


P(//)=P(//IF)xP(F)+P(//IB)xP(F)=(0.50)(0.90)+(0.70)(0. 10)=0.52 . (38) 

Applying Bayes’ theorem, we can update our prior probabilities P(F) and P(B) after we have 
observed 10 consecutive heads as follows: 

Define event of 10 consecutive heads as H\q. Thus we obtain: 

P(R\H l f^ojgjfW (a7 10 XftlOJ 0 763 ™, 

Similarly, we obtain for: 


P(F\Hio)=l-P(B I //i 0 )=0.237 . (40) 

We observe that the experiment has resulted in an increase of the probability that the coin is biased 
and a corresponding decrease that it is “honest.” As mentioned before, the real problem lies in the 
assignment of the prior probabilities. 

NOTE: Many objections to Bayes’ theorem are actually attacks on Bayes’ postulate (also called the 
Principle of Indifference or Principle of Insufficient Reason), that says if we have no knowledge of the 
prior probabilities, we may assume them to be equally probable. In our case we would set 
P(B)=P(F)= 0.5, which is, of course, a very dubious assumption. 

AUXILIARY QUESTIONS: (Solutions are left as a challenge to the reader.) 

(1) A military operation consists of two independent phases. Phase A has a 10-percent risk and 
phase B, a 20-percent risk. (Risk is defined as the probability of failure.) What is the probability of 
mission failure? Answer: 0.28 . 


(2) In a straight poker hand, what is the probability of getting a full house? Answer: • 

(3) Two prizes are awarded in a lottery consisting of 200 tickets. If a person buys two tickets, 
what is the probability of winning the first prize or the second prize, but not both (exclusive “OR”)? 
Answer: 0.0198995 . 

(4) A student recognizes five different questions that may be asked on a quiz. However, he has 
time to study only one of them that he selects randomly. Suppose the probability that he passes the test if 
“his” question appears is 0.90, but the probability that he passes the test if it does not appear is only 0.30. 
The test contains only one question and it is one of the five. 

(a) What are the chances that he will pass the test? Answer: 0.42 . 
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(b) If the student passed the test, what is the probability that “his” question was asked on the 
test? Answer: 0.4286 . 

(5) A man has two pennies — one “honest” and one two-headed. A penny is chosen at random, 
tossed, and observed to come up heads. What is the probability that the other side is also a head? 

Answer: 2/3 . 


D. Probability Distributions 

In order to define probability distributions precisely, we must first introduce the following 
auxiliary concepts. 

1. Set Function 

We are all familiar with the concept of a function from elementary algebra. However, a precise 
definition of a function is seldom given. Within the framework of set theory, however, it is possible to 
introduce a generalization of the concept of a function and identify some rather broad classification of 
functions. 

An ordered pair consists of two elements, say a and b, in which one of them, say a, is designated 
as the first element and the other the second element. 

The Cartesian product AxB of two sets A and B is the set of all ordered pairs (a, b ) where aeA 

and bsB. For instance, if A=( a, b, c) and B=(l, 2, 3, 4) then the Cartesian product AxB is illustrated by 
the following area. Only the enclosed area F represents a function as defined below in figure 5. 



A relation R from a set A into a set B is a subset of AxB. 


A function F from a set A into a set B is a subset of AxB such that each ae A appears only once 
as the first element of the subset (see fig. 6). 
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Function F 



The domain is the set of all elements of A which are related to at least one element of B. 

The range is the set of all elements of B which are related to at least one element of A. 

EXAMPLE: Supermarket: 

A=Products (domain): x 
B = Price (range): y => y-y (*). 

A set function is a function where the elements of the domain are sets and the elements of the 
range are numbers. 

2. Random Variable 

It is often desirable to assign numbers to the nonnumerical outcomes of a sample space. This 
assignment leads to the concept of a random variable. A random variable X is a set function which 
assigns to each outcome EsS a real number X(E)=x. 

The domain of this set function is the sample space S and its range is the set of real numbers. In 
general, a random variable has some specified physical, geometrical, or other significance. It is important 
to observe the difference between the random variable itself and the value it assumes for a typical outcome. 

It has become common practice to use capital letters for a random variable and small letters for the 
numerical value that the random variable assumes after an outcome has occurred. Informally speaking, we 
may say that a random variable is the name we give an outcome before it has happened. It may be that the 
outcomes of the sample space are themselves real numbers, such as when throwing a die. In such a case, 
the random variable is the identity function X(E)=E. 

Note that, strictly speaking, when we are throwing a die the outcomes are not numerical, but are 
the “dot patterns” on the top face of the die. It is, of course, quite natural, but really not necessary, to 
associate the face values with the corresponding real number. 

3. Probability Function and Cumulative Distribution Function 

EXAMPLE: A coin is tossed three times. The sample space then consists of eight equally probable out- 
comes: S= { HHH. . . 777} . Let X be the random variable that counts the number of heads in each outcome. 
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Thus, X has range {0, 1, 2, 3}. Figure 7 lists the value x that the random variable X assumes for each 
outcome and the probability associated with each value x. 



Note that different outcomes may lead to the same value of x. 

A probability function is a function that assigns to each value of the random variable X the 
probability of obtaining this value: 

f(x)=P(X=x) 0<f(x)<l . (41) 

For example, /(l)=P(X=l)=3/8. 

The mathematical function for/(x) is (see fig. 8). 



A cumulative distribution function (c.d.f.) is a function that assigns to each value of the random 
variable X the probability of obtaining at most this value. It is sometimes called “or less” cumulative 
distribution: 


P(X<x)='Lf x (t)=F x (x) . (42) 

t<x 

Occasionally the complementary cumulative distribution function (c.c.d.f.) is also called the “or 
more” c.d.f. is used (i.e., nuclear power industry). One is easily obtained from the other. 
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Because of its appearance, the step function is also called staircase function (see fig. 9). Note that 
the value at an integer is obtained from the higher step, thus the value at 1 is 4/8 and not 1/8. As we 
proceed from left to right (i.e., going “upstairs”) the distribution function either remains the same or 
increases, taking values between 0 and 1 . Mathematicians call this a monotonically increasing function. 


m 


4/8 A 


*r 


1/8 H 


1 ► x 

0 12 3 

Figure 9. — Cumulative distribution function. 


In general f{x) is a probability function if: 


o</M<i 

(43) 

3£m 

ii 

(44) 


4. Continuous Random Variable and Probability Density Function 

Examples of continuous random variable and probability density function are temperature, age, and 
voltage. These ideas can be extended to the case where the random variable X may assume a continuous 
set of values. By analogy we define the cumulative distribution function for a continuous random variable 
by: 

P(X<x)= J f x (t)dt = F x (x) . (45) 


The rate of change of P(X<x) with respect to the random variable (probability/interval) is called 
the probability density function (p.d.f.): 

• ( 46 , 

Note that the probability density function is not a probability. The probability density function has the 
following properties: 

f(x)> 0 (47) 


\f(x)dx = 1 . (48) 
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Also note that the probability of a null event (impossible event) is zero. The converse is not neces- 
sarily true for a continuous random variable, i.e., if X is a continuous random variable, then an event 
having probability zero can be a possible event. For example, the probability of the possible event that a 
person is exactly 30 years old is zero. 

Alternate definitions: 


P(x<X<x+dx)-f(x) dx 

(49) 

b 

P(a<X<b)=jf(x)dx=F(b)-F(a) . 

(50) 


a 


E. Distribution (Population) Parameters 

One of the purposes of statistics is to express the relevant information contained in the mass of data 
by means of a relatively few salient features that characterize the distribution of a random variable. These 
numbers are called distribution or population parameters. We generally distinguish between a known 
parameter and an estimate thereof, based on experimental data, by placing a hat symbol ( A ) above the 

estimate. We usually designate the population parameters by Greek letters. Thus, fl denotes an estimate of 
the population parameter fi. In the following, the definitions of the population parameters will be given 
for continuous distributions. The definitions for the discrete distribution are obtained by simply replacing 
integration by appropriate summation. 

1. Measures of Central Tendency (Location Parameter) 

a. Arithmetic Mean (Mean, Average, Expectation, Expected Value). The mean is the most 
often used measure of central tendency. It is defined as: 


fi= jxf(x)dx . (51) 


The definition of the mean is mathematically analogous to the definition of the center of mass in dynamics. 
That is the reason the mean is sometimes referred to as the first moment of a distribution. 

b. Median (Introduced in 1883 by Francis Galton). The median m is the value that divides 
the total distribution into two equal halves, i.e., 


m 


F(m)= \f(x )dx - 1/2 

— oo 


(52) 


c. Mode (Thucydides, 400 B.C., Athenian Historian). This is also called the most probable 
value, from the French “mode,” meaning “fashion.” It is given by the maximum of the distribution, 
i.e., the value of x for which: 


f£W =0 

ax 


(53) 
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If there is only such maximum, the distribution is called unimodal; if there are several, it is called 
multimodal. 

In a symmetrical unimodal distribution, the mean, median, and mode coincide (see fig. 10). The 
median and the mode are useful measure for asymmetric (skewed) distributions. The median is commonly 
used to define the location for the distribution of family incomes, whereas the mode is used by physicians 
to specify the duration of a disease. 



Figure 10. — Location of mean, median, and mode. 

It is useful to remember that the mean, median, and mode of a unimodal distribution occur in 
reverse alphabetical order for a distribution that is skewed to the right, as in figure 10, or in alphabetical 
order for a distribution that is skewed to the left. What gives the mean the great importance in statistical 
theory is its mathematical tractability and its sampling properties. 


2. Measures of Dispersion 


The degree to which the data tend to spread about the mean value is called the dispersion or 
variability. Various measures of this dispersion are given below: 

(a) Variance (Introduced by R.A. Fisher in 1918): 

o 2 =l(x-v) 2 f(x)dx=\x 2 f(x)dx-Li 2 =E(x 2 )-{E(x)} 2 . (54) 

(b) Standard deviation (Introduced by K. Pearson in 1890): The standard deviation is simply 
the positive square root of the variance, denoted by < 7 . 

(c) Mean deviation: 


o d =\\x-ii\f(x)dx . (55) 

For a normal distribution, the mean deviation =4/5 cr. 
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(d) Coefficient of variation (c.o.v.): The c.o.v. gives the standard deviation relative to the 
mean fi as: 

c.o ,v.=(T/fi . (56) 

It is a nondimensional number and is often expressed as a percentage. Note that it is independent of the 
units used. 


3. Quantiles 

The quantile of order p is the value i; p such that P(X<% p )=p. 

For p= 0.5 Median 

p=0.25 Quartile 
/?=0.10 Decile 
p=0.01 Percentile. 

The y'th quantile is obtained by solving for 

j/n= jf(x)dx . (57) 

— oo 

Quantiles are sometimes used as measures of dispersion: 

(a) Interquartile range <2=£o.75-£o.25 

(b) Semi-interquartile range <22=0.5 x (<^0.75-^0.25) 

(c) Interdecile range <2 10=^0.90-^0.10 

For a normal distribution, the semi-interquartile range =2/3 o. 


4. Higher Moments 

Other important population parameters referred to as “higher moments” are given below: 

(a) Skewness: 

= where =j( x—flftf(x)dx . (58) 

A distribution has a positive skewness (is positively skewed) if its long tail is on the right and a negative 
skewness (is negatively skewed) if its long tail is on the left. 

(b) Kurtosis: 

cc 4 =p 4 /<T 4 ,'wheren 4 =l(x-n) 4 f(x)dx . (59) 
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Kurtosis measures the degree of “peakedness” of a distribution, usually in comparison with a normal 
distribution, which has the kurtosis value of 3. 

(c) Moments of kth order: 

Moments about the origin (raw moments): 


H' k =E(x k )=jx k f(x)dx . (60) 

Moments about the mean (central moments): 

H k = E(x-ji) k =j(x-/J.) k f(x)dx . (61) 

F. Chebyshev’s Theorem 

If a probability distribution has the mean p and the standard deviation a, the probability of 

obtaining a value that deviates from the mean by at least k standard deviations is at most *//c 2 (see 
fig. 11). Expressed in mathematical form: 

P(l X-p \>k o)> l/k 2 (upper bound) (62) 

P(\X-p I < k a)>\-l/k 2 (lower bound) . (63) 



Figure 11. — Chebyshev’s theorem. 


Chebyshev’s theorem is an example of a distribution-free method; i.e., it applies to any distribution 
with the only condition that its mean and standard deviation exist. This is its strength but also its 
weakness, for it turns out that the upper and lower bounds are too conservative to be of much practical 
value (see table 3 for a comparison with the k-factors derived from the normal distribution). 


TABLE 3. — Normal distribution compared with Chebyshev’s theorem. 


Percentage 

Normal JT-factor 

Chebyshev ff-factor 

90.0 

1.65 

3.16 

95.0 

1.96 

4.47 

99.0 

2.58 

10.0 

99.9 

3.29 

31.6 
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The so-called Camp-Meidel correction for symmetrical distributions is only a slight improvement. It 
replaces the upper bound by 1/(2. 25x£ 2 ). 


G. Special Discrete Probability Functions 


1. Binomial Distribution 

Assumptions (Bernoulli trials) of binomial distribution are: 

• There are only two possible outcomes for each trial (“success” and “failure”) 

• The probability p of a success is constant for each trial. 

• There are n independent trials. 

The random variable X denotes the number of successes in n trials. The probability function is then 
given by: 

f(x) = P(X = x)=fyp x (\-p) n - x forjt=0,l...n (64) 

Mean fl=np 

Variance <y 2 =npq where q=\-p 
Skewness 0 : 3 = %=£^ 

sjnpq 


Kurtosis <% 4 = 3 + l r — ~ . 

4 n pq 

(Note that any of the subsequent BASIC programs can be readily edited to run on any computer.) 

Binomial probabilities can be calculated using the following BASIC program, which is based on 
the recursion formula: 


f(x)=£x^xf(x-l) 


(65) 


10 

15 

20 

25 

30 

35 

40 

45 

50 

55 


"CUMBINOMIAL" 

CLEAR: INPUT "N="; N, "P="; P,"X="; X 

Q=1-P:F=Q A N:B=F:S=0 

IF X=0 GOTO 55 

FOR 1=1 TO X 

E=P*(N-I+ 1 )/Q/I 

F=F*E:S=S+F 

NEXT I 

CF=S+B 

PRINT CF: PRINT F 
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EXAMPLE: N= 6, P=0.30, jc=3 


CF= 0.9295299996, F=1 .852199999E-01 . 

Simulation of Bernoulli trials can be performed using the following BASIC program: 

100: "BINOMIAL" 

105: INPUT "N=";N, "P=";P 

110: S=0 

1 15: FOR 1=1 TO N 

120: U=RND .5 

125: IF U<P LET X=1 GOTO 135 

130: X=0 

135: S=S+X 

140: NEXT I 

145: PRINTS: GOTO 110 


2. Poisson Distribution 

The Poisson distribution is expressed as: 

for x=0, 1,2,3 ...oo (66) 

Mean: [£=p 


Variance: 


The Poisson distribution is a limiting form of the binomial distribution when n-+° o, p-+0, and 
np = jU remain constant. The Poisson distribution provides a good approximation to the binomial 
distribution when n>20 and p<0.05. 

The Poisson distribution has many applications that have no direct connection with the binomial 
distribution. It is sometimes called the distribution of rare events; i.e., “number of yearly dog bites in New 
York City.” It might be of historical interest to know that its first application by Ladislaus von Bortkiewicz 
( Das Gesetz der Kleinen Zahlert, Teubner, Leipzig, 1898) concerned the number of Prussian cavalrymen 
killed by horse kicks. 

The following BASIC program is based on the recursion formula: 

f(x)=%f(x- 1) . (67) 


10: "POISSON" 

15: INPUT "M=";M, "X=";X 

20: P=EXP-M:S=1:Q=1 
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25: IF X=0 LET F=P GOTO 50 

30: FOR 1=1 TO X 

35: Q=Q*M/I: S=S+Q 

40: NEXT I 

45: F=P*S:P=P*Q 

50: PRINT F:PRINT P 

60: END 


EXAMPLE: fi= 2.4 x=4 

F=9.041314097xl0->, P=/(jc)=1.2540848986x10- 1 . 


3. Hypergeometric Distribution 

This distribution is used to solve problems of sampling inspection without replacement, but it has 
many other applications (i.e., Florida State Lottery, quality control, etc.). 

Its probability function is given by: 


f(x)= 


N-a 

n—x 


for x=0, 1, 2, 3 ...n 


( 68 ) 


where N=lot size, n=sample size, a=number of defects in lot, and jc=number of defects (successes) 
in sample. 


Mean: /u = nxj^ 


Variance: a 1 = 


nxax(N-a)x(W-n) 

W 2 x(N-1) 


(69) 


The following BASIC program calculates the hypergeometric probability function and its 
associated cumulative distribution by calculating the logarithms (to avoid possible overflow) of the three 
binomial coefficients of the distribution and subsequent multiplcation and division, which appears in line 
350 as summation and difference of the logarithms. 


300 

305 

310 


"HYPERGEOMETRIC" 
INPUT "N=";N, "Nl, "A=";A 
INPUT "X=";X1 


315: NC=N:KC=Nl:GOSUB 400 

320: CD=C 


325: CH=0: FOR X=0 TO XI 

330: NC=A:KC=X:GOSUB 400 

335: C1=C 

340: NC=N-A:KC=Nl-X:GOSUB 400 

345: C2=C 

350: HX-EXP (C1+C2-CD) 

355: CH=CH+HX 

360: NEXT X:BEEP3 
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365: PRINT CH: PRINT HX 

370: END 

400: C=0: ;M=NC-KC 

405: IF KC=0 RETURN 

410: FOR 1=1 TO KC 

415: C=C+LN((M+I)/I) 

420: NEXT I : RETURN 


4. Negative Binomial (Pascal) Distribution 

The negative binomial, or Pascal, distribution finds application in tests that are terminated after a 
predetermined number of successes has occurred. Therefore, it is also called the binomial waiting time 
distribution. Its application always requires a sequential-type situation. The random variable X in this case 
is the number of trials necessary for k successes to occur. (The last trial is always a success!) 

The probability function is readily derived by observing that we have a binomial distribution up to 
the last trial, which is then followed by a success. As in the binomial distribution, the probability 
p of a success is again assumed to be constant for each trial. 

Therefore: 


k for x=k, £+1, k+2 ... (70) 

where r=Number of trials 

/c=Number of successes 
/?=Probability of success. 

Also: 


Mean: p=— Variance: <y 2 = ^ . (71) 

P p z 

The special case £=1 gives the geometric distribution. The geometric distribution finds application in 
reliability analysis, for instance, if we investigate the performance of a switch that is repeatedly turned on 
and off, or a mattress to withstand repeated pounding in a “torture test.” 

PROBLEM: What is the probability of a mattress surviving 500 thumps if the probability of failure as the 
result of one thump is p= 1/200? 

SOLUTION: The probability of surviving 500 thumps is P(X>500) and is called the reliability of the 
mattress. The reliability function R(x) is related to the cumulative distribution as: 

R(r)=l-F(r) . (72) 

Since the geometric probability function is given as f (x)=p( \-p) , we obtain for the reliability of 

the mattress: 
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DO 


( 73 ) 


R(x)= 'Lf(x)=(l-p) 500 =( 0.995) 500 = 0.0816 . 
jr=501 

It is sometimes said that the geometric distribution has no “memory” because of the relation: 

P(X>x 0 +x \X>xo)=P(X>x) . (74) 

In other words, the reliability of a product having a geometric (failure) distribution is independent of its 
past history, i.e., the product does not age or wear out. These failures are called random failures. 

Sometimes the transformation z=x-k is made where z is the number of failures. Then the 
probability function is given by: 


f(z) = {^ +Z z l ^P k ^q z for z=0, 1,2,3... (75) 

and these are the successive terms in p k {\-q)~ z , a binomial expression with a negative index. 

When a program of the binomial distribution is available, the negative probabilities can be obtained 
by the simple identity: 


f p (x\k,p)=^f b (k\x,p) (76) 

where the subscript p denotes the Pascal distribution and the subscript b the binomial distribution. 

PROBLEM: Given a 70-percent probability that a job applicant is made an offer, what is the probability of 
needing exactly 12 sequential interviews to obtain eight new employees? 


SOLUTION: 

*=12, &=8, p=0.7 (/i=l 1 .42, <t=2.2) 

(77) 


f p ( 1218, 0.7=( 1 ? 1 )(0.7) 8 (0.3) 4 =0.1541 . 

(78) 


A more interesting complementary relationship exists between the cumulative Pascal and the 
cumulative binomial distribution. This is obtained as follows: The event that more than x sequential trials 
are required to have k successes is identical to the event that x trials resulted in less than k successes. 
Expressed in mathematical terms we have: 

P p (X>x I k,p)=P b (K<k I x,p) . (79) 

The cumulative Pascal distribution can then be obtained from the cumulative binomial by the 
relationship: 


F p {x\k,p)=\-Fb(k-\ \x,p) . (80) 

PROBLEM: Given a 70-percent probability that a job applicant is made an offer, what is the probability 
that, at most, 15 sequential interviews are required to obtain eight new employees? 

SOLUTION: F p (\5 1 8, 0J)=1-F b (7 1 15, 0.7)=0.95=95 percent . (81) 
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H. Special Continuous Distributions 


1. Normal Distribution 

The normal distribution is the most important continuous distribution for the following reasons: 

• Many random variables that appear in connection with practical experiments and observations are 
normally distributed. This is a consequence of the so-called Central Limit Theorem or Normal 
Convergence Theorem to be discussed later. 

• Other variables are approximately normally distributed. 

• Sometimes a variable is not normally distributed, but can be transformed into a normally 
distributed variable. 

• Certain, more complicated distributions can be approximated by the normal distribution. 

The normal distribution was discovered by De Moivre (1733). It was known to Laplace no later 
than 1774, but is usually attributed to Carl F. Gauss. He published it in 1809 in connection with the theory 
of errors of physical measurements. In France it is sometimes called the Laplace distribution. 
Mathematicians believe the normal distribution to be a physical law, while physicists believe it to be a 
mathematical law. 

The normal distribution is defined by the equation: 


g(y) = 2 v y , 0 < y < <» . f(x) = 

fy 


i 


0\2tc<j 




for — °o < x < oo 


(82) 


Mean: fi-ji Variance: <fi=cP- 

Skewness: 0:3=0 Kurtosis: 04=3. 

The cumulative distribution function of this density function is an integral that cannot be evaluated 
by elementary methods. 

The existing tables are given for the so-called standard normal distribution by introducing the 
standardized random variable. This is obtained by the following transformation: 


4 a ' 


(83) 


In educational testing, it is known as the standard score or the “z-score.” It used to be called 
“normal deviate” until someone perceived that this is a contradiction in terms (oxymoron), in view of the 
fact that deviates are abnormal. 


Every random variable can be standardized by the above transformation. The standardized random 
variable always has the mean jl=0 and the standard deviation cx=l . 
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The cumulative distribution of every symmetrical distribution satisfies the identity: 

F(-x)=\-F(x) . 


(84) 



The error function is defined as: 


0(x)=~r= \e dt . 

V7T 0 

It is related to the normal cumulative distribution by: 



■ 

r v Yl 

1+0 

X 


v \ 2 J_ 


where 


F(x)= ^ Je u 2/2 du . 

\ ~~.lt — oo 


(85) 


( 86 ) 


(87) 


We are often interested in the probability of a normal random variable to fall below (above) or 
between k standard deviations from the mean. The corresponding limits are called one-sided or two-sided 
tolerance limits. Figures 12 and 13 are given as illustrations. 


Two-Sided Limits 
Pr ( n -Ko<X<fi +Ko)=A 



FIGURE 12. — Normal distribution areas: two-sided tolerance limits. 


Stated, there is a 95.46-percent probability that a normal random variable will fall between 
the ±2 a limits. 
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One-Sided Limits 
Pr(X< g + Ko) = F 



FIGURE 13. — Normal distribution areas: one-sided tolerance limits. 


Stated, there is a 97.73-percent probability that a normal random variable will fall below (above) 
the +2(7 (-2(7) limit. 


The relationship between the two areas A and F is: 

A=2F-1 or F=(A + l)/2 . 


A normal random variable is sometimes denoted as: 

X=i V(/z, cr 2 )=Gauss (jU, a 1 ) . 


( 88 ) 

(89) 


The standard normal variable or standard score is thus: 

Z=N(0, l)=Gauss (0, 1) . (90) 


Table 4 gives the normal scores (“A'-factors”) associated with different levels of probability. 


TABLE 4. — Normal K-factors. 


One-Sided 


Two-Sided 


Percent 

Ki 

Percent 

Kz 

99.90 

3.0902 

99.90 

3.2905 

99.87 

3.0000 

99.73 

3.0000 

99.00 

2.3263 

99.00 

2.5758 

97.73 

2.0000 

95.46 

2.0000 

95.00 

1.6448 

95.00 

1.9600 

90.00 

1.2815 

90.00 

1.6449 

85.00 

1.0364 

85.00 

1.4395 

84.13 

1.0000 

80.00 

1.2816 

80.00 

0.8416 

75.00 

1.1503 

75.00 

0.6744 

70.00 

1.0364 

70.00 

0.5244 

68.27 

1.0000 

65.00 

0.3853 

65.00 

0.9346 

60.00 

0.2533 

60.00 

0.8416 

55.00 

0.1257 

55.00 

0.7554 

50.00 

0.0000 

50.00 

0.6745 
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As was already mentioned, the normal p.d.f. cannot be integrated in closed form to obtain the 
c.d.f.. One could, of course, use numerical integration, such as Simpson’s Rule or the Gaussian 
quadrature method, but it is more expedient to calculate the c.d.f. by using power series expansions, 
continued fraction expansions, or some rational (Chebyshev) approximation. An excellent source of 
reference is Handbook of Mathematical Functions by Milton Abramowitz and Irene A. Stegun (eds.). 

The following approximation for the cumulative normal distribution is due to C. Hastings, Jr.: 


F(x)= 1 - 


e ~* 2/2 

-\12n 


\ a n y n 

n = 1 


(91) 


1 

1 + 0.231 64 19jc 


0 < x < °° 


(92) 


where 

ai=0.3193815 

a 2 =-0.3565638 

03=1.781478 Error=<7E-7 

a 4 =-l. 821256 

a 5 =1.330274. 

EXAMPLE: x=2.0 F(*)=0.97725 


Sometimes it is required to work with the so-called inverse normal distribution, where the area F is 
given and the associated X-factor (normal score) has to be determined. A useful rational approximation for 
this case is given by equation (93) (equation 26.2.23 in the Handbook of Mathematical Functions) as 
follows: 


We define Q(k p )=p where Q=\-F and 0 <£><0.5 . 


Then 


K p =t- 


Cft+Cit + Cjt , v ^ „ — - 

0 1 2 -j + £(p), t = V~2 in p 


l + d x t + d 2 t + dft 


(93) 


and Ie(p)l<4.5xl0- 4 , 
where 

c 0 =2.515517 z/j = 1 .432788 

Cj=0.802853 rf 2 =0.189269 

c 3 =0.01 0328 ^3=0.001308. 

A more accurate algorithm is given by the following BASIC program. It calculates the inverse 
normal distribution based on the continued fraction formulas also found in the Handbook of 
Mathematical Functions. Equation (94) in this reference (equation 26.2.14 in the Handbook ) is used for 
jc>2, equation (95) for jc<2 (equation 26.2.15 in the Handbook, and the Newton-Raphson method using 
initial values obtained from equation (98). 
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(94) 


x> 0 

*>0 (95) 

x =/ ^ — ^—~+£(t), t= / l«-i-, and \e(t) I<3xl0 - ^ (96) 

p l+£>,/+£> 2 t z V / 2 

a 0 =2.30753 £>] =0.99229 

a 1=0.27061 £>2=0.0448 1. 

NORMIN (0.5<P<1) 

DEFDBL A-Z 

INPUT"P=";P:PI=3. 141592653589793# 

REM Equation 26.2.22 

Q=l-P:T=SQR(-2*LOG(Q)) 

A0=2.30753:A1=.27061:B1=.99229:B2=.0481 
NU=A0+A 1 *T:DE= 1 +B 1 *T+B2*T*T 
X=T-NU/DE 


L0: Z=l/SQR(2*PI)*EXP(-X*X/2) 
IF X>2 GOTO LI 


REM Equation 26.2.15 

V=-25-13*X*X 

FOR N=ll TOO STEP- 1 

U=(2*N+ 1 )+(- 1 ) A (N+ 1 )*)N+ 1 )*X*X/V 

V=U:NEXT N 

F=.5-Z*X/V 

W=Q-F:GOTO L2 

REM Equation 26.2.14 
LI: V=X+30 
FOR N=29 TO 1 STEP-1 
U=X+N/V 
V=U:NEXT N 
F=Z/V:W=Q-F:GOTO L2 

REM Newton-Raphson Method 
L2: L=L+1 
R=X:X=X-W/Z 
E=ABS(R-X) 

IF E>.0001 GOTO L0 

PRINT USING "##.####" ;X 
END 
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The normal distribution is often used for random variables that assume only positive values such as 
age, height, etc. This can be done as long as the probability of the random variable being smaller than 
zero, i.e., P(X< 0), is negligible. This is the case when the coefficient of variation is less than 0.3. 

The normal and other continuous distributions are often used for discrete random variables. This is 
admissible if their numerical values are large enough that the associated histogram can be reasonably 
approximated by a continuous probability density function. Sometimes a so-called continuity correction 
is suggested which entails subtracting one-half from the lower limit of the cumulative distribution and 
adding one-half to its upper limit. 

2. Uniform Distribution 


The uniform distribution is defined as f(x)= T -—, a<x<b . The mean of the uniform distribution 

J ' ' b-a 

is fi = — and the variance <7 2 Y1 ' • The uniform p.d.f. is depicted in figure 14. 


A 

m 

i 

b-a 


a 


b 


+- 

x 


Figure 14. — Uniform p.d.f. 


3. Log-Normal Distribution 

Many statisticians believe that the log-normal distribution is as fundamental as the normal 
distribution itself. In fact, by the central limit theorem, it can be shown that the distribution of the product 
of n independent positive random variables approaches a log-normal distribution, just as the sum of n 
independent random variables approaches a normal distribution. It has been applied in a wide variety of 
fields including social sciences, economics, and especially in reliability engineering and life testing. The 
log-normal distribution is simply obtained by taking the natural logarithm of the original data and treating 
these transformed data as a normal distribution. 


In short, Y- in X is normally distributed with log mean fly and log standard deviation ay. Since 
we are really concerned with the random variable X itself, we have to determine the probability density 
of X. By the methods shown later, it can be shown that it is: 


f(x) = 


1 

xo Y y2Jt 



/ \ 2 

£nx~{l Y | 

v ) 


, x >0 


(97) 


It can also be shown that f(0) = f'(0)=0. 
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The mode of the distribution occurs at: 


and the median at: 


The mean is: 


The variance is: 


* - e VY-°Y 

x Mode ~ e 

X —p^Y 

^Median ~~ e 
H x =e M Y +(l/2)a Y 


(98) 

(99) 


( 100 ) 


a 2j e 2n r +a 2 y 



( 101 ) 


The distribution has many different shapes for different parameters. It is positively skewed, the 
degree of skewness increasing with increasing Oy. 

Some authors define the log-normal distribution in terms of the Briggsian (Henry Briggs, 1556- 
1630) or common logarithm rather than on the Naperian (John Napier, 1550-1616) or natural logarithm. 

The log mean fly and the log standard deviation ay are nondimensional pure numbers. 


4. Beta Distribution 


This distribution is a useful model for random variables that are limited to a finite interval. It is 
frequently used in Bayesian analysis and finds application in determining distribution-free tolerance limits. 

The beta distribution is usually defined over the interval (0, 1) but can be readily generalized to 
cover an arbitrary interval (jcq, * 1 ). This leads to the probability density function: 


f(x)= 


1 T (oc+p) 

(x x -x Q )T(a)T(p) 


*-*o f ^ 

x l~ x o) 



*-*o 

x l~ x o) 


( 102 ) 


The standard form of the beta distribution can be obtained by transforming to a beta random 
variable over the interval (0, 1) using the standardized z-value: 


z = — — , where 0<z<l . 
x l~ x 0 


(103) 


The beta distribution has found wide application in engineering and risk analysis because of its 

diverse range of possible shapes (see fig. 15 for different values of a and /?.) The mean and variance of 
the standardized beta distribution is: 


/*= 


a 

a+fi 


and <r 2 = 


a/3 

(a+f3) 2 (a+f}+l) 


(104) 
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Figure 15. — Examples of standardized beta distribution. 


The cumulative beta distribution is known as the incomplete beta function. An interesting and 
useful relationship exists between the binomial and beta distribution. If X is a binomial random variable 
with parameters p and n, then: 


P(X<x) = 


Tfn+lJ 


T(n-x)Y(x+\) 


1 ~P 

jt n ~ x ~ ] a-t) x 

0 


dt 


(105) 


5. Gamma Distribution 

The gamma distribution is another two-parameter distribution that is suitable for fitting a wide 
variety of statistical data. It has the following probability density function: 

r a-l p -x//5 

f(x)=- — for x>0 and a, B > 0 . (106) 

J r(a)P a H 

A typical graph of the gamma p.d.f. is shown in figure 16. 
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The parameter (3 is a scale parameter that only changes the horizontal scale. The parameter a is 
called the index or shape parameter because changes in a result in changes in the shape of the graph of the 
p.d.f. The quantity T(a) represents the gamma function defined by: 

oo 

r(cc)=jx a - l e- x dx . (107) 

0 

Integration by parts shows that T(a)=(a-1) r(a-l). If a is a positive integer then T(a)=(a-1)! 
The gamma distribution has many applications in inventory control, queuing theory, and reliability studies. 
In reliability work it is often found to be in competition for application with the Weibull distribution, which 
is discussed in the following section. 


The moment generating function of the gamma distribution is given by: 

M( t)=E[e tx Ml-pt)~ a 


(108) 


from which the first two central moments can be derived as: 

Mean: H=af3 

Variance: a 2= a /j2 . (109) 

This gamma distribution includes the exponential distribution as a special case when 0=1. The 
exponential distribution is commonly used in the study of lifetime, queuing, and reliability studies. Its 

single parameter /3 is the mean of the distribution and its inverse is called the failure rate or arrival rate. In 
reliability studies the mean is often referred to as the mean time between failures (MTBF). 

Another special case of the gamma distribution is obtained by setting oc=v/2 and /3=2. It is called 

the chi-square (% 2 ) distribution, which has many statistical applications as discussed subsequently in more 
detail. 

When the parameter a is restricted to integers, the gamma distribution is also called the Erlangian 
distribution. In this case, the cumulative distribution can be obtained by integration by parts. Otherwise, 
the gamma density cannot be integrated in closed form and the cumulative probability distribution is then 

referred to as the incomplete gamma function. The Erlangian represents the waiting time distribution for a 
independent events, each of which has an exponential distribution with mean (3. 

6. Weibull Distribution 

This distribution was known to statisticians as the Fisher-Tippett Type III asymptotic distribution 
for minimum values or the third asymptotic distribution of smallest extreme values. It was used by the 
Swedish engineer Waloddi Weibull (1887-1979) in the analysis of the breaking strength of materials 
(1951). It has gained quite some popularity in reliability engineering due to its mathematical tractability and 
simple failure rate function. It is extensively used in life testing where it is competing with the gamma and 
log-normal distributions in search for the true nature of the random variable under investigation. The 
Weibull distribution can be easily memorized using its cumulative distribution function, which is: 

F(x)= l-e-vtP =l-e~( x/T l) P (110) 


where r\ is called characteristic life. 
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It is easily seen that for x=rj the cumulative distribution is always 1-1/^=0.6321 regardless of the 
value of the parameter /3. The first form of the cumulative distribution is somewhat easier to manipulate 

mathematically. If the parameter a is known, the characteristic life can be calculated as r\=a - y/ A The 
probability density is obtained by simple differentiation of the cumulative distribution. 

There are two special cases that give rise to distributions of particular importance: 

CASE 1: /J=l 

This is the exponential distribution, which is also a special case of the gamma distribution. It has many 
important applications in reliability and life testing. 

CASE 2: P=2 

This is the Raleigh distribution. It describes the distribution of the magnitude of winds if the north-south 
winds and the east-west winds are independent and have identical normal distributions with zero mean. It 
is also the distribution for the circular error probability (CEP). Another application arises in random 
vibration analysis where it can be shown that the envelope of a narrow-band random process has a 
Rayleigh distribution. 

7. Extreme Value (Gumbel) Distribution 

This distribution is seldom used in reliability engineering and in life and failure data analysis. It is, 
however, used for some types of “largest observations” such as flood heights, extreme wind velocities, 
runway roughness, etc. It is more precisely called the largest extreme value distribution. Sometimes it is 
also called die Gumbel distribution, after Emil J. Gumbel (1891-1967), who pioneered its use. It is briefly 
presented here for the sake of completeness. 

Its cumulative distribution function is: 

F(x)=exp{-exp (-(x-A)/<5) } , -oo <X <°° (111) 

Mean: p=X+y^, where y=0. 57722 (Euler’s constant) (112) 

Variance: <7 2 =7t 2 S 2 /6 . (113) 


I. Joint Distribution Functions 


1. Introduction 

In many engineering problems we simultaneously observe more than one random variable. For 
instance, in the cantilever beam shown in figure 17, we encounter five random variables: the length L, the 

force P, the Young’s modulus E, the deflection S, and the moment of inertia I z . In general, all of them 
will, of course, have different distributions as indicated in figure 17. 
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FIGURE 17. — Cantilever beam. 


DEFINITION: A A>dimensional random vector X=[X i, X 2 ...Xk\ is a set function that assigns to each 
outcome EeS real numbers X, (E)=Xj (i=l,2...k). 

2. Discrete Case 

The bivariate event is expressed as: 

p(x l ,x 2 )=P{(X l =xi) n (X 2 =x 2 )) . (114) 

This is called the (bivariate) joint probability function. Extensions to higher dimensional random vectors 
are self-evident. The bivariate probability function has the obvious properties: 

p(x i,*2)>0 and £ 'Lp(x l ,x 2 )=l . (115) 

* 1*2 

The cumulative bivariate distribution is defined as: 

F(x i,x 2 )= I L P(u v u 2 ) . (116) 

M l -*1 « 2-*2 


A natural question arises as to the distribution of the random variables X\ and X 2 themselves. 
These are called marginal probability functions and are defmed as: 

P\(x\)=lP(x v x 2 ) (117) 

*2 


and 


P2(x2)='LP(x 1 ,x 2 ) . (118) 

*l 

NOTE: The terminology stems from the fact that if the probabilities associated with the different events are 
arranged in tabular form, the marginal probabilities appear in the margin of the table. 

The conditional probability function is defined as: 

PI (X 1 ,x 2 >= P< p X 2 (xj ) > ■ < 119 > 
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3. Continuous Case 


The extension to continuous random variable is easily achieved by replacing summation by 
integration (and usually by replacing p with/). 

DEFINITION: The c.d.f. is defined as: 


x \ x 2 

F(XpX 2 )=P(X l <*], X 2 <x 2 ) = J jf(u l ,u 2 )du i du 2 


and the bivariate probability density is obtained by partial differentiation as: 


f(x v x 2 )= 


d 2 F(x\,x 2 ) 
dx^ dx 2 


Similarly, the conditional probability density is defined as: 


( 120 ) 


( 121 ) 


M ■ (122) 

4. Independent Random Variables and Bayes’ Rule 

The notion of independence, which was introduced in connection with intersections of events, can 
be analogously defined for random vectors. Informally speaking, when the outcome of one random 
variable does not affect the outcome of the other, we state that the random variables are independent. 

DEFINITION: If/(jq, xj) is the joint probability density of the random variables X\ and A/, these 
random variables are independent if and only if: 

/(*l,*2)=/l(*l)/2(*2) • (123) 

Extensions to the multivariate case are obvious. 

Sometimes the concept of independence is defined in terms of the cumulative distributions. In this 
case we have: 


F(x\, xi)=F\(xi) F 2 (x 2 ) . (124) 

Bayes’ Rule. In recent years there has been a growing interest in statistical inference, which looks 
upon parameters (mean, standard deviations, etc.) as random variables. This approach is known as 
Bayesian estimation. It provides a formal mechanism for incorporating our degree of personal judgment 
about an unknown parameter into the estimation theory before experimental data are available. As these 
data become available, the personal judgment is constantly updated to reflect the increase in knowledge 
about a certain situation. 

DEFINITION 1 : The probability distribution h 0 (p) of the parameter p, which expresses our belief 
about the possible value p before a sample is observed, is called the prior distribution of p. We reiterate 
that h 0 (p) makes use of the additional subjective information about the value p prior to taking a sample. 
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DEFINITION 2: The updated probability distribution of the parameter p, which expresses our 
increase in knowledge, is called posterior distribution. This distribution is calculated from the 
relationship that exists between the joint and conditional probability distributions as: 

f(x,p)=g\(x\p) h 0 (p)=g 2 (p\x) h 2 (x) . (125) 

The posterior distribution is then given by: 


g 2 (p\x) = 


g^}P)hg(Pl 

tyx) 


(126) 


where /i 2 (x)=marginal distribution. Here g\{x \ p) is the sampling distribution of x given the parameter p. 

It is readily recognized that the above formula is Bayes’ theorem as applied to the continuous case. 
Once the posterior distribution is obtained, it can be used to establish confidence intervals. These 
confidence intervals are called Bayesian confidence intervals. 

It should be mentioned that if the prior information is quite inconsistent with the information gained 
from the sample, the Bayesian estimation yields results that are inferior to other estimation techniques; 
e.g., the maximum likelihood estimation. The Bayesian estimation technique has become a favorite in 
many different statistical fields. 

As an application, we discuss the estimation of the unknown parameter p of a binomial 
distribution: 

g l (x\p)^p x (l-p) n - x ( 127 ) 


«=N umber of trials 
jt=Number of successes 
/^Probability of success. 

We assume the prior distribution of p to be the beta distribution: 

*o (p)= rwrw p(a ~'’ v-pf™ • < 128 > 

We first calculate the marginal distribution of x for the joint density fix, p) by integrating over p: 

1 1 

h 2 (x) = jf(x,p)dp=jg l (x\p)h 0 (p)dp . (129) 

0 0 

This yields: 



T(a+p) 

r (a)T(P) 


l 


jp( x+a 0(1— p/ n X+ P ^ dp 
0 


(130) 


or: 
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h ( r \-( n '\ r («+ft) T(a+x)r(n-x+P) 

^ {X) UJ r(a)T(p) r (n+a+P) 

The posterior distribution is, therefore: 


/ i \ — r (n+cx+P) 
g 2 (p x)- r(a+x)r(N _ x+ p) 


pi x+a-l ) _ p )(n-x+P-l) 


(131) 


(132) 


We have established the following theorem: If x is a binomial distribution and the prior distribution 
of p is a beta distribution with the parameters a and /?, then the posterior distribution of g 2 (p I x) is also a 
beta distribution with the new parameters x+a and n-x+fi. 

The expected value (mean) of the posterior distribution is: 


1 


P = jP 82 (P u ) d P = 
0 


X+GC 

n+a+b 


(133) 


EXAMPLE 1: Uniform Prior Distribution: If the prior distribution is uniform, i.e., ho(p)= 1, then a=(5-l 
and the mean of the posterior distribution is given by: 


P = 


x+1 

n+2 


(134) 


This is known as the Laplace law of succession. It is an estimate of the probability of a future event given 
that x successes have been observed in n trials. Sometimes this law is stated as follows: The more often an 
event has been known to happen, the more probable it is that it will happen again. 

EXAMPLE 2: Tests With No Failures: Assuming a uniform prior distribution, we calculate the posterior 
distribution for a situation where n tests were performed without a failure. We obtain the posterior 
distribution g 2 (p I n)=(n + 1 ) p n (see fig. 18). 
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p = 0.50 (Mean) 





FIGURE 18. — Posterior distribution with no failures. 

EXAMPLE 3: Two Tests and One Failure: The posterior distribution is gjip I x)=6p (1 -p) (see fig. 19). 
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EXAMPLE 4: Bayesian Confidence Limits: In reliability studies we are usually interested in one-sided 
confidence limits for the reliability p of a system. Given the posterior distribution g 2 (p I x) for the reliability 
p, we obtain the lower Bayesian (1-a) confidence limit as: 

1 

\g 2 (p\x)dp=\-a . (135) 


The lower confidence limit is obtained by solving this equation for the lower integration limit pp 
(see fig. 20). 



FIGURE 20. — Lower confidence limit. 


The upper Bayesian (1-a) confidence limit is similarly obtained from the equation: 

Pu 

jg 2 (p\x)dp=l-a (136) 

0 

If there are no test failures the posterior distribution is: 

g 2 (p\x)=(n+l) p n . (137) 

Inserting this posterior distribution into the above equation for the lower confidence limit for the 

reliability (probability of success), we can solve it for the (1-a) lower confidence limit in closed form and 
obtain: 


p (n+ l) =a ( 138 ) 

This lower confidence limit is often referred to as the demonstrated reliability when it is applied to a 
system which has undergone n tests without a failure. 

It is sometimes necessary to find the number of tests required to demonstrate a desired reliability 

for a (1 - a) confidence level. This can be easily obtained by solving the above equation for n, which in 
this case is: 


47 



(139) 


in a 

n=en "L 


-1 


For instance, if it is required to demonstrate a 99-percent reliability at the 50-percent confidence level, it is 
necessary to run 68 tests without a failure. 

If no confidence level is specified, the 50-percent confidence level is generally assumed. This is 
especially the case for high reliability systems. 

The non-Bayesian confidence limit for the binomial parameter p for the above case is given by: 

P L n =a (140) 


It is somewhat more conservative than the Bayesian limit because it requires one more test to demonstrate 
the same level of reliability. 


J. Mathematical Expectation 

In many problems of statistics, we are interested not only in the expected value of a random 
variable X, but in the expected value of a function Y=g(X). 

DEFINITION: Let X=(X], X 2 ...X„) be a continuous random vector having a joint probability density 
f {x)—f{x j , X 2 -- X n ) and let u(X) be a function of X. The mathematical expectation is then defined as: 


E[u{X)]= ju(x)f(x)dx . (141) 


For a discrete random vector the integral is replaced by a corresponding summation. Sometimes 
£[.] is called the “expected value operator.” 

The vector notation introduced here simplifies the expression for the multiple integral appearing in 
the definition of the mathematical expectation. A further consideration of this general case is beyond the 
scope of the present text. 

EXAMPLE: Chuck- A-Luck: The concept of a mathematical expectation arose in connection with the games 
of chance. It serves as an estimate of the average gain (or loss) one will have in a long series of trials. In 
this context the mathematical expectation is also called the expected profit. It is defined by the sum of all 
the products of the possible gains (losses) and the probabilities that these gains (losses) will occur. 
Expressed mathematically, it is given as: 


E=la iPi . (142) 

It is important to keep in mind that the coefficients a, in the above equation are positive when they 
represent gains and negative when they represent losses. 

In betting odds, the expected profit is set equal to zero (£= 0). For example, if the odds of favoring 
an event are 3:1, the probability of the event happening is 3/4 and the expected profit is: 

£=$1 (3/4)-$3 (l/4)=0 . (143) 
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In Scame’s New Complete Guide to Gambling, John Scarae puts it this way: “When you make a 
bet at less than the correct odds, which you always do in any organized gambling operation, you are 
paying the operator the percentage charge for the privilege of making a bet. Your chances of winning has 
what is called a ‘minus expectation.’ When you use a system you make a series of bets, each with a minus 
expectation. There is no way of adding minuses to get a plus.” 

As an illustration for calculating the expected profit, we take the game of Chuck-A-Luck. In this 
game the player bets on the toss of three dice. He is paid whatever he bets for each die that shows his 
number. 

The probability of any number showing on a die is p= 1/6. Let the probability of one number 
showing on three dice bepj, of two numbers p 2 , and of three numbers p^. Then: 



pi =3 pq 2 = 3 (1/6) (5/6) (5/6)=75/216 

(144) 


P2=3p 2 <7=3 (1/6) (1/6) (5/6)= 15/2 16 

(145) 


p3=p 3 =(l/6) (1/6) (1/6)=1/216 . 

(146) 

The expected profit is, therefore: 



£=($ 1 ) p 1 +($2) p 2 +($3) p3— ($ 1 ) P LOSS 

(147) 

where 

P LO ss = 1 ~(P 1 + P2+P3)= 1 25/2 1 6 . 

(148) 

Thus, 

c 1x75+2x15+3x1 1w 125 17 

216 1X 216 — 216° 

(149) 

or: 

£=-7.87 cents/dollar. 

(150) 


The house percentage is, therefore, 7.87 percent. 

As may be noticed, we have already encountered a few special cases of mathematical expectation, 
such as the mean, variance, and higher moments of a random variable. Beyond these, the following simple 
functions are of particular interest. 

1. Covariance 

The covariance measures the degree of association between two random variables X\ and X 2 . It is 
defined as follows: 


Cov (Xu X 2 )=<Ji 2 =E[(X l -m) (X 2 -P2)] . 
Performing the multiplication we obtain an alternate form of the covariance as: 

Cov (Xi, X 2 )=£[Xi xX 2 ]-£[Xi] E[X 2 ]=E[X 1 xX 2 l-pi p 2 . 


(151) 


(152) 
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The covariance can be made nondimensional by dividing by the standard deviation of the two 
random variables. Thus, we get the so-called correlation coefficient: 


Cov(X v X 2 ) <7j2 
slV(X = ' 


(153) 


It can be shown that if X\ and Xj are independent, then p= 0. This is true for both discrete and 

continuous random variables. The converse of this statement is not necessarily true, i.e., we can have p= 0 
without the random variables being independent. In this case the random variables are called uncorrelated. 

It can also be shown that the value of p will be on the interval (-1, +1), that is: 


-l<p <+l 


(154) 


2. Linear Combinations 


The linear combination of the elements of a random vector X=(Xi, X 2 -..X n ) is defined as 

u(X)—Y—ao+a\ X+...+a/ iX n . (155) 

The mean and variance of this linear combination is: 


and: 


E(Y) = fl Y = a 0 + 'La i E(X i )=a 0 +J t a,p, 
/= 1 1=1 


n n 


V(Y)=Oy = 'Lafof + X Sfljfl iOn 
1=1 1=1 j = 1 




(156) 


(157) 


If the variables are independent, the expression for the variance of Y is simplified because the covariance 
terms <7y are zero. In this case the variance reduces to: 


<r?=Ia?<7? . (1 58) 

i=l 

Notice that all these relationships are independent of the distribution of the random variables 
involved. 


K. Functions of Random Variables 

In many engineering applications we meet with situations in which we have to determine the 
probability density of a function of random variables. This is especially true in the theory of statistical 
inference, because all quantities used to estimate population parameters or to make decisions about a 
population are functions of the random observation appearing in a sample. 
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For example, suppose the circular cross section of a wire is of interest. The relationship between 
the cross section of the wire and its radius is given by A=nR 2 . Considering R a random variable with a 
certain distribution, the area A is also a random variable. Our task is to find the probability density of A if 
the distribution of R is known. 

DEFINITION: Let B be an event in the range space R x , i.e., BeR x , and C be an event in range space R y , 
i.e., C € R y such that: 

B={xeR x :h(x)eR y } . (159) 

Then B and C are equivalent events — that means they occur simultaneously. Figure 21 illustrates a 
function of a random variable. 



FIGURE 21 . — A function of a random variable. 


Equivalent events have equal probability: 


P{A)=P{B)=P{C) . (160) 

1. Discrete Random Variables 

The p.d.f. of a function of a discrete r.v. is determined by the equivalent events method. 

EXAMPLE: In the earlier coin-tossing problem, the random variable X can assume the four values 
x=0, 1, 2, 3 with respective probabilities 1/8, 3/8, 3/8, 1/8. Let the new variable be Y=2X- 1. Then the 
equivalent events are given by y=- 1, 1, 3, 5 with probabilities p Y (-l)=l/8, p Y (l)=3/8, p Y (3)=3/8, 
p Y (5)=l/8. Notice that equivalent events have equal probabilities. 

2. Continuous Random Variables 

The p.d.f. of a function of a continuous r.v. can be determined by any of the following methods: 
Method I: Transformation of Variable Technique (“Jacobian” Method). 

Method II: Cumulative Distribution Function Technique (Equivalent Event Method). 

Method HI: Moment Generating Function Technique. 

APPLICATIONS: 

Method I: 

(a) Univariate case. Let y=h(x) be either a decreasing or an increasing function. Then: 

f(x) dx=g(y) dy, (161) 
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dx 

dy 


or solving for g(y): 

g(y)=f(*) | 

Note that the absolute value sign is needed because of the condition that g(y) > 0. 
EXAMPLE 1: Let: f(x)=e~ x for 0<x<°o . 

The new variable is given as y=x 2 / 4. 


(162) 


(163) 


Thus: 


Therefore: 


^-=x/2^>~ = 2/x . 
dx dy 


g(y)=e 


-X 


=—e~ x 

x 


(164) 


(165) 


In terms of the new variable y: g(y) = — —e ~ '^ y ,0 < y < °o . 


(166) 


EXAMPLE 2: Random sine wave: y=A sin jc (see figs. 22 and 23). The amplitude A is considered to be 
constant and the argument jc is a random variable with a uniform distribution: 


f(x)=± with-f <x<f 


(167) 
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(fy 

dx 


= Acosx, 


dx _ 1 

dy A cos x 


>0 ,g(y)= 


1 1 
;r cosx 


(168) 


Expressed in terms of the new variable y : 

cos( x) = y 1 - sin 2 (x) = yj 1 -( y/A ) 2 (169) 

and finally: g(y)= r ~ * . (170 

nA^l-fy/A) 1 



FIGURE 23. — Probability density of random sine wave. 


EXAMPLE 3: Random number generator (probability integral transformation): A random variable X is 
given that has a uniform distribution over the interval (0, 1). Find a new random variable Y that has a 
desired probability distribution g(y). 

According to method I, we have the relationship: 

f(x) dx=g(y) dy . (171) 

Since f(x)-l for 0<r<l we obtain from this: dx=g(y) dy 


and after integration, 

y 

x= jg(u)du = G(y) (172) 

— oo 

where, of course, G(y) is the cumulative distribution of y. Note that the cumulative distribution F(x) has a 
uniform distribution over the interval (0, 1) independent of/(x). 
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The desired random variable Y is then obtained by solving the above equation for the inverse 
cumulative function, such that 

y=G _1 (x), where G ~ 1 (x) is often referred to as the percentile function. (173) 

Figure 24 shows the relationship between the two random variables in graphical form. 



EXAMPLE: Determine the transformation law that generates an exponential distribution, such that: 

g(y)=ae~ ce y 0 <y <°° (174) 


G(y)= joce au du= 1-e 
0 

x=G(y)= l-e-W 


y=-±£n(l-X) . 


Random numbers can then be generated from: 


1 


Bivariate case: 


y;=-a tnx i 


f(x, y) dxdy=g(u, v) dudv 


(175) 

(176) 

(177) 

(178) 


As is known from advanced calculus, the differential elements dxdy and dudv are related through 
the Jacobian determinant as: 
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dudv=J(u, v)dudv where 7(«,v)=Jacobian determinant. 


(179) 


dxdy = 


dx 

dx 

du 

dv 

dy 

dy 

du 

dv 


The new joint probability density function is, then: 

g(u, v)=f(x, y) I /(«, v) I . (180) 

EXAMPLE 1: If the random variables X\ and X 2 have a standard normal distribution N(0, 1), the ratio 
X\/X 2 has a Cauchy distribution. 


Let y\=x\ix 2 and y2=*2 => *l=yi/y2 and X2=y2- 
The partial derivatives are: 


dx 


1 _ 


<fy\ } ’ 2 ’ dy 2 7V ^ 


dx^ 


=y\ 


— =0, -=^=\ 


dy 2 


(181) 


and the Jacobian is then J=y2- 


The joint probability density function is: 


g(yi>y 2 ) 


2n 



hxi W\=^' iyl,Uy 'U 


(182) 


The marginal distribution is: 


8(y ' )= l 

which gives: g(y x )= . 

n(Uyf) 

EXAMPLE 2: Box-Muller method: An important use of the Jacobian method for finding the joint 
probability density of a function of two or more functions of a random vector is the Box-Muller algorithm. 
Since its introduction in 1958, it has become one of the most popular method for generating normal 
random variables. It results in high accuracy and compares favorably in speed with other methods. The 
method actually generates a pair of standard normal variables from a pair of uniform random variables. 

Consider the two independent standard normal random variables x and y with densities: 


V 2 7 


(183) 

(184) 


f(x)= 


' 


.1*2 1 -iv 2 

2 and /fy)= 1 e. 


din 


(185) 
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Let us introduce the following polar coordinates: 

x=r cos 0 and y=r sin 0 . (186) 

The joint density of the polar coordinates can be obtained from the joint probability density of the Cartesian 
coordinates using the Jacobian method as: 


g(r, <}>)=f(x, y) 171 

where 171 is the absolute value of the Jacobian determinant: 


J= 


dx 

dr 

dx 

d<p 


COS( 0) 

-rsin(0) 

dy_ 

dy_ 


sin(0) 

rcos(0) 

dr 

d<p 





= rcos z ('0)+rsin 2 ('0)=r 


Because of the independence of the Cartesian normal random variables we have: 

r2 


f(x,y)=f(x)f(y)=^e 1 2 (x2+y2) 


=J-e~2 

2n e 


and finally using the value of the Jacobian determinant: 


(187) 


(188) 


(189) 


g(r,<t>)=^ e 


(190) 


Since the new joint distribution is a product of a function of r alone and a constant, the two polar 
random variables r and 0 are independent and have the following distribution: 


g l (r)=re 2 , O<0<°° (191) 

and: 

g 2 (0)=2^ , O<0<27r . (192) 

The random variable r has a Rayleigh distribution that is a special case of the Weibull distribution 
with a=l/2 and /3=2. It can be simulated using the probability integral transformation. The random 
variable 0 is a uniform distribution over the interval (0, 2k). Therefore, we can now generate the pair of 
independent normal random variables x and y using two uniform random numbers u\ and «2 by the 
following algorithm: 


x = ^'-2 in («j ) cos(2ku 2 ) and y = -J-2 £n(u x ) sin( 2k u 2 ) (193) 


Method II: 

IfXi,X 2 ...X„ are continuous random variables with a given joint probability density, the 
probability density of y=h{x\, X 2 - ■ -x n ) is obtained by first determining the cumulative distribution: 
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G(y)=P(Y<y)=P [h(x \,X2,... x n )<y ] 


(194) 


and then differentiating (if necessary) to obtain the probability density of y: 

g ( y ) = ~, j^T (195) 

EXAMPLE 1: Z=X + Y. This is one of the most important examples of the function of two random 
variables X and Y. To determine the cumulative distribution F(z) we observe the region of the xy plane in 
which x+y<z. This is the half plane left to the line defined by x+y=z as shown in figure 25. 



FIGURE 25. — Sum of two random variables. 


We now integrate over suitable horizontal strips to get: 


oo z-y 

F (z)= I J f xy (x,y)dxdy 


(196) 


— oo — oo 


Assuming independence of the random variables, we have fxyix , y)=f x (x)fy(y). Introducing this in 
the above equation, we obtain: 


\z-y 


F (z) = J 1 1 fr(x)dx 


-oo — oo 


■f v (y)dy= \ F x (z-y)fjy)dy 


(197) 


Differentiating with respect to z, we get the p.d.f. of z as: 


f(z)= j f x (z-y)f y (y)dy= j f x (x)f (z-x)dx . (198) 

— oo — oo 

Note that the second integral on the right-hand side is obtained if we integrate over vertical strips in the left 
half plane. 
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RESULT: The p.d.f. of the sum of two independent random variables is the convolution integral of their 
respective probability densities. 

EXAMPLE 2: Z=X-Y. Another important example is the probability of the difference between two random 
variables. Again we find the cumulative distribution F(z ) by integrating over the region in which x-y < z 
using horizontal strips as indicated in figure 26. Following similar steps as in the previous example, we 
obtain the cumulative distribution of Z as: 


F(z)= ! 


z+y 

i f xy (x,y)dxdy . 


(199) 


y 



FIGURE 26. — Difference of two random variables. 


Assuming independence of the random variables as above yields: 


F(z)= 


z+y 

j f x (x)dx 


/v 


(y)dy= J F x (z+y)f y (y)dy 


( 200 ) 


Differentiating with respect to z, we get the p.d.f. of z as: 


f(z)= \ f x (z+y)f y (y)dy= j f x (x)f y (z+x)dx . (201) 


RESULT: The p.d.f. of the difference between two independent random variables is the correlation 
integral of their respective probability densities. 

APPLICATION: Probabilistic failure concept. One way a system can fail is when the requirements 
imposed upon the system exceed its capability. This situation can occur in mechanical, thermal, electrical, 
or any other kind of system. The probabilistic failure concept considers the requirement and the capability 
of the system as random variables with assumed probability distributions. Let C represent the capability 
(strength) having a distribution fdx) and R the requirements (stress) of the system having a distribution 
fR(y). Then system failure occurs if R>C or, in other words, if the difference U-C-R is negative. The 
difference U is called the interference random variable (see fig. 27). 
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Interference Area 


FIGURE 27. — Interference random variable. 


The cumulative distribution of U is given as: 


F(u) = 


oo 


(u+y 


OO 



J f c (x)dx 


\fR(y)dy= J F c (u+y)f R (y)dy 

-oo 


( 202 ) 


Therefore, the probability of failure, also called the unreliability of the system, is: 


F(0)= j 


y 

j f c (x)dx 

-oo 


f R (y)dy= J F c (y)f R (y)dy . 


(203) 


In general, the integration has to be done numerically. However, if the two random variables are 
normally distributed, we know that the difference of the two random variables is again normally 
distributed because of the reproductive property of the normal distribution. 

Introducing the normal score of the difference u, we have: 


^ / 2 1 


(204) 


The probability of failure is then given by setting u=0. By design, the probability of failure will be 
very small, such that the normal score will turn out to be negative. To use the normal score as a measure of 

reliability, it is therefore customary to introduce its negative value, which is called the safety index (3. It is 
numerically identical with the one-sided AT-factor and the above equation then becomes: 


ff=~ p 2 = • (205) 

V a C+°R 

The corresponding reliability can then be obtained from the normal A'-factor table or using a 
program that calculates the cumulative normal distribution. A good reference book on this subject is 
Reliability in Engineering Design by Kapur and Lamberson. 
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Method III: 


This method is of particular importance for the case where we encounter linear combinations of 
independent random variables. The method uses the so-called moment generating function. 

DEFINITION: Given a random variable X the moment generating function of its distribution is given by: 


M(t)= p( X: ) for discrete X (206) 

/=! 


and 

oo 

M{t)= \e tx f(x)dx for continuous X . (207) 

— oo 

The moment generating function derives its name from the following moment generating property: 


and 


4HMH 1 = y X n e tx , p ( x ) for discrete X 
dt n i=\ 1 1 


(208) 


dfM(t} 

dt n 


\x n e tx f(x)dx 

— oo 


Taking the derivatives and setting t = 0 gives: 


for continuous X . 


and 


d n M(0) 

dt n 


= I xfp(x i )=p n 
i = 1 


for discrete X 


(209) 


( 210 ) 


— ^(0) _ j x n^ x ^ ( j x for continuous X . (211) 

dt n _oo 

Readers who are familiar with the Laplace transform will quickly notice that the moment generating 
function for a continuous random variable X is identical with the definition of the (two-sided or bilateral) 
Laplace transform if one substitutes t=-s. In fact, if one defines the moment generating function as the 
bilateral Laplace transformation of the probability density function, it is possible to use existing Laplace 
transform tables to find the moment generating function for a particular p.d.f.. 

In some advanced textbooks the variable t in the moment generating function is replaced by “i t," 

where i= \/— 1 • It is then called the characteristic function of the distribution. This measure is sometimes 
necessary because certain distributions do not have a moment generating function. 

Method HI makes use of the following theorem: 

If Xi, X 2 ...X„ are independent random variables and T=Xi+X 2 +...X m , then 


60 



(212) 


n 

My(t)=TlM x Jt) . 

1=1 ' 

EXAMPLE: The normal distribution has the moment generating function: 


M n (t)=e 


pt+j<y 2 t 2 


(213) 


It is readily seen that the sum of k normal random variables is again a normal distribution with 

mean and variance =£<r? . This is called the reproductive property of the normal 

distribution. Only a few distributions have this property. The Poisson distribution is another one that has 
this property. 


L. Central Limit Theorem (Normal Convergence Theorem) 

If X\, X 2 ,...X n are independent random variables with arbitrary distributions such that 


mean: 

E(xi)=pi t 

(214) 

and variance: 

V(xi)=o? , 

(215) 

then: 

i/W 

(216) 

has an approximate standard normal distribution jV( 0, 1) for large n. 
importance and ubiquity of the normal distribution. 

This is the basic reason for the 


APPLICATION: Normal random number generator: 


12 

Gauss(0,l)=N(0,l)=X£/.(0,l)-6 . (217) 

i=l 

Summing up 12 random numbers, which are uniform over the interval (0,1), gives an easy method 
and, for practical purposes, accurate standard normal random variable. It is truncated at ±6 < 7 . The 
probability of exceeding these limits is less than 2x1 0 -9 . 

M. Simulation (Monte Carlo Methods) 

Simulation techniques are used under the following circumstances: 

• The system cannot be analyzed by using direct and formal analytical methods. 

• Analytical methods are complex, time-consuming, and costly. 

• Direct experimentation cannot be performed. 

• Analytical solutions are beyond the mathematical capability and background of the 
investigator. 
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EXAMPLE: Buffon’s needle: This famous problem consists of throwing fine needles of length L on a 
board with parallel grid lines separated by a constant distance a. To avoid intersections of a needle with 
more than one line, we let a>L (see fig. 28). 


A 


a 

i 

Figure 28. — Buffon’s needle. 



We define the position of a needle in terms of the distance x of its midpoint to the nearest line, and 
its direction with respect to the grid line by the angle 6. 

For an intersection to occur, we have the condition: 

x<L/2 cos 0 (218) 


where the range of the angle is Q<B<Jt!2 and the range of the distance 0<x<aJ2. 

Both random variables x and 6 are assumed to be uniform over their respective interval. 
Geometrically, we can observe that this condition is given by the shaded area under the cosine curve in 
figure 29. 

The rectangle represents all possible pairs of (x,0) and the shaded area all pairs (x,ff) for which an 
intersection occurs. Therefore, probability of an intersection is given by the ratio of these two areas. The 

area under the cosine curve is A\=L/2 and the area of the rectangle is An— a tcIA. Taking this ratio, we 
obtain the probability as: 


/> =V A 2=5S • (219) 



62 



The following BASIC program simulates Buffon’s needle experiment using equal needle length 
and grid distance (a=L= 1). For this condition the probability of an intersection is obtained from the above 
formula as p=0.6362. 

Using a run size of 10,000 for the simulation, we can calculate the maximum error of estimate, 
which is given by: 


E= 


z a/2 Aj 


\Pjl~P) 


( 220 ) 


Using an 80-percent confidence level for which z a /2=1.28, we have a maximum error of 
approximately 0.006. The corresponding confidence interval is then: 


0.6302</?<0.6422 . 


( 221 ) 


BUFFON’S NEEDLE 

PI=4*ATN(1) 

RANDOMIZE 123 


L1:S=0: FOR 1=1 TO 10000 

X=RND 

Y=COS(PI/2*RND) 

IF X<Y THEN X=l:GOTO L2 
X=0 

L2:S=S+X:NEXT I 

BEEPrBEEP 
PRINT S:GOTO LI 
END 


The computer simulation gave the following 12 results: /?=0.6315, 0.6325, 0.6355, 0.6352, 
0.6382, 0.6382, 0.6353, 0.6230, 0.6364, 0.6357, 0.6364, 0.6453, 0.6362. We observe that 2 out of 12 
runs fall outside the 80 percent confidence interval, which is what one would expect. 

An actual experiment was supposed to have been carried out by an Italian mathematician named 
Lazzarini in 1901, who made 3,408 needle tosses and counted 2,169 intersections. This gave a value for 

71=3.142461964, which is only wrong in the third decimal place. One can easily recognize that it is highly 

unlikely to achieve such an accuracy with only 3,408 tosses. It is more likely that Lazzarini knew that n is 
approximately Tin . If we substitute this number in the above expression for the probability, we get: 

P _ 2x7 _ 14x155 _ 2170 

22 22x155 “3410 ‘ ^ 222 ) 

Of course, this might have been too obvious, so the numbers were slightly altered. What would have been 
the result if he had thrown one more needle? 
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III. STATISTICS 


A. Estimation Theory 

The concepts of probability theory presented in the preceding chapters begin with certain axioms (a 
self-evident fact that itself is not proved, but which is the basis of all other proofs) and derive probability 
laws for compound events. This is a deductive process. In statistics we are dealing with the inverse 
process, which is an inductive one: given certain observations, we want to draw conclusions concerning 
the underlying population. The study of such inferences is called statistical inference. Probability looks 
at the population and predicts a sample. Statistics looks at the sample and predicts the population. 

1. Random Sample 

In order that conclusions drawn from statistical inference are valid, it is necessary to choose 
samples that are representative of the population. The study of sampling methods and concomitant 
statistical analysis is called the design of experiments. 

One can easily see that it is often difficult to get a random sample. For example, if we are after a 
random sample of people, it is not good enough to stand on a street comer and select every 10th person 
who passes. Many textbooks present the following classical example of a faulty sampling method. In 
1936, a magazine called Literary Digest sent its readers and all telephone subscribers — approximately 
10 million people — a questionnaire asking how they intended to vote at the forthcoming presidential 
election. 

There were 2,300,000 replies, from which it was confidently predicted that the Kansas Republican 
Governor Alfred M. Landon would be elected. As it turned out, he lost every state but Maine and Vermont 
to the incumbent president Franklin D. Roosevelt. 

This erroneous prediction happened because readers of this literary magazine and people who had 
telephones did not, in 1936, represent a fair, i.e., random sample of the American voters. It was especially 
risky to overlook the enormous proportion of nonreplies. 

DEFINITION: If X\, X 2 . . . X n are independent and identically distributed (iid) random variables, they are 
said to constitute a random sample of size n. The joint distribution of this set of random variables is: 

g(xi ,X 2 -. .x n )=f{x\)f(x 2 ).. .f{x n ) (223) 

where f(x) is the population distribution. 

The term population is used to describe the sample space (universe) from which the sample is 
drawn. Some statisticians use the term parent population. 

2. Statistic 

A statistic is a random variable that is a function of the random variables that constitute the random 
sample. The probability distribution of a statistic is called sampling distribution. 
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3. Estimator 


An estimator is a statistic that is used to make a statistical inference concerning a population 
parameter ft 


Q=Q(X v X 2 ...X n ) 


(224) 


B. Point Estimation 

The value of an estimator provides a point estimate of a population paramater 0. 

1. Properties of Estimators 

There are often many possible functions capable of providing estimators. In the past the choice of 
an estimator was sometimes governed by mathematical expediency. With the arrival of high-speed 
computational means, this aspect is no longer as dominant as it used to be. For example, in on-line quality 
control the range was commonly used to estimate the dispersion of data rather than the standard deviation. 

We list several desirable characteristics for good estimators as follows: 

(a) Unbiasedness: 


E(©)=6 


(225) 


(b) Consistency (asymptotically unbiased): 


lim p(|©-0<e} = l or lim E(Q)=G 
n—>oo LI J rt— »o o 


(226) 


(c) Efficiency: If 


Var(0 1 )<Viflr(0 2 ) , (227) 

^ -A. 

0j is more efficient than 0 2 if both are unbiased. 

(d) Sufficiency. An estimator is sufficient if it extracts all the information in a random sample 
relevant to the estimation of the population parameter ft Mathematically expressed as the Fisher-Newman 

A A 

factorization theorem, 0 = 0 (x\, X 2 -.-x n ) is sufficient if the joint p.d.f. of X\,X 2 ...X n can be factorized 
as: 


A 

/(*!, X2-..x n \0)=g[Qxi,X2 ...x n I 0]h(x i,X 2 ... x n ) (228) 

where g depends on x\, x 2 -..x n only through © and h is entirely independent of ft 

EXAMPLE 1: Let X\, X 2 -..X n be normally distributed. By setting /l= 0 \ and a 2 - 62 and 6 ={ 0 \, 62), we 
have: 
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f(x v x 2 ... x n \0)= exp 

[\ 2n& 2 J 


■ 2 ^|, (x r e > >2 


(229) 


so that 


'L(x j -G l ) 2 = '£(x j -x) 2 +n(x-O l ) 7 
j = 1 7 j = 1 7 




- 4 - 


2tt0 2 J 2y 2 7=1 7 2« 2 


(230) 


(231) 


It follows that X, X ( X • -X) 2 is sufficient for 0=(0i, # 2 ). 

I 7=1 7 J 


EXAMPLE 2: Consider the uniform p.d.f.: 


/(x)=^forO<jc < 0 


=0 otherwise 


(232) 


where 0 is to be estimated on the basis of n independent observations. We want to prove that the statistic 
0=max (jcj, X2 -.. x n ) is a sufficient estimator of the parameter 0. 


SOLUTION: It can be shown that the cumulative distribution function of 0 is 


V'HfF 


for Oot<0 


(233) 


and, therefore, the p.d.f. of 0 is /§(•*)= for O<x<0. Therefore, since the joint p.d.f. of Xj, 


X 2 -..X„ is , which may be factored as 


( 4 ) =[ j — 7-7 > or f( x u x 2 ■•■x n \Q)=g[®( x \’ x 2---Xn\Q)\h(xi,X2...x n ) (234) 

\0/ v 0 J nx n ~ l 

as was to be shown. It is seen that the right-hand side of this expression has the desired product form such 
that the first term is a p.d.f. of the statistic © and the second term does not depend on the parameter 6. 

(e) Mean square error (MSE). The MSE of an estimator is defined as the expected value of the 
square of the deviation of the estimator from the parameter © being estimated. It can be shown that the 
MSE is equal to the variance of the estimator plus the square of its bias: 
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MSE=E(Q-d) 2 = E[G-E(e)] 2 +[E(e)-6] 2 . (235) 

A biased estimator is often preferred over an unbiased estimator if its MSE is smaller. It is often possible 
to trade off bias for smaller MSE. 

A A 

PROBLEM: Suppose O i and 02 are two different estimators of the population parameter 0. In addition, 

we assume that E(Q i)=0, £(02)=O.9 0, Var (0 j)=3 and Var (02>=2. Which estimator has the smaller 
MSE? 

SOLUTION: MSE ( 0 ] )=Var ( 0 i )+(Bias) 2 =3+0=3 (236) 

MSE ( 0 2 )=Var ( © 2 )+(Bias) 2 =2+0.0 1 0 2 . (237) 

We notice that as long as 1 0 1 < 10 the estimator ©2 has a smaller MSE and is, therefore, “better” 

A 

than the estimator 0 1 , in spite of the fact that it is a biased estimator (see fig. 30). 



FIGURE 30. — Sampling distribution of biased and unbiased estimator. 


EXAMPLE 1: Variance of normal distribution. A commonly used estimator for the variance of a 
distribution is the statistic: 


s2 . ZXj-X ) 2 
n-1 


(238) 


The popularity of this estimator stems from the fact that it is an unbiased estimator for any distribution. 
However, it is important to realize that the corresponding sample standard deviation s, which is the 

positive square root of the variance a 2 , is a biased estimator of the population standard deviation a. 
Another frequently used estimator of the variance of a population is the statistic: 


o2_»*,-XJ 2 


(239) 


If the sample is taken from a normal distribution, the above estimator is, in fact, the maximum likelihood 
estimator, but it is not unbiased. 
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The MSE of the two estimators, when calculated for a normal distribution with variance a 2 , is: 

MSEf5 2 j=^- and MSE(sl)= (2n ~l - <7 - . (240) 

Tl — 1 D fi ^ 

It can be readily verified that the MSE of the biased estimator is smaller than that of the unbiased estimator 
for all sample sizes n. For a sample size of n= 5, the MSE of the unbiased estimator is 39 percent higher 
than that of the biased one. 

EXAMPLE 2: Binomial parameter p. As a particular interesting example we compare two estimators 
available for the parameter p of the binomial distribution. The unbiased estimator of this parameter is: 



(241) 


with x being the number of successes and N the number of trials. This is sometimes called the “natural” 
estimator and is, coincidentally, also the maximum likelihood estimator. The other biased estimator is the 
Bayesian estimator obtained by using a uniform prior distribution and is given by: 


A _ X+\ 

p b~ N+2 


(242) 


It can be shown that the variance of the unbiased estimator is Var ( p u )=^jj- and the variance of the 
Bayesian estimator Var ( p b )~ • Clearly the variance of the biased estimator is smaller than the 

unbiased one. To obtain the MSE we need to determine the bias of the Bayesian estimator. Its value is 
given by: 


Bias= 


l-2p 

W+2 ' 


(243) 


Figure 31 shows the absolute value of the bias of the estimator as a function of the parameter p. 

It is seen that the bias is zero for p=0.5 and is a maximum for p=±l. Theoretically, the MSE of the 
biased estimator is smaller than that of the unbiased one for: 

pq> 4(M+l) • (244) 
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FIGURE 3 1 . — Estimator bias as a function of parameter. 


For /V— > oo, the critical value for /?=0.854, implying that for greater values the unbiased estimator 
should be preferred because of its smaller MSE. 

However, as the following BASIC Monte Carlo program demonstrates, the Bayesian estimator 
will always be better for the practical case regardless of the value of p. The discrepancy between the 
theoretical case and the real case stems from the fact that the theoretical variance of both estimators 
becomes very small for high values of the parameter p. As a consequence, the bias term becomes dominant 
and tilts the inequality in favor of the unbiased estimator. 

'BINOMIAL ESTIMATOR (4/2/92) 

BEGIN: 

INPUT"N=",N 
INPUT"P=",P 
RANDOMIZE 1 2 3 

'HERE STARTS THE MONTE CARLO RUN WITH N=100 
L2:EB=0:EC=0: FOR R=1 TO 100 

THE NEXT 5 LINES ARE THE BERNOULLI TRIALS 

S=0:FOR B=1 TO N 

U=RND 

IF U<=P THEN X=l:GOTO LI 

x=o 

L1:S=S+X:NEXT B 

'CALCULATION OF MSE 
PC=X/N :PB=(X+ 1 )/(N+2) 

EC=EC+(PC-P) A 2:EB=EB+(PB-P) A 2 

NEXT R: BEEPrBEEP 

PRINT EC,EB 
GOTO L2 
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Sample run: N=5; p- 0.90 

MSE ( Pu ) =53.16, 52.20, 53.16, 54.12, 50.28 


MSE ( Pb ) =40.28, 39.69, 40.28, 40.87, 38.52 

It is seen that the MSE of the Bayesian binomial estimator is always smaller than that of the 
classical estimator. 


2. Methods of Point Estimation 

There are three methods of point estimation that will be discussed: 

• Method of (matching) Moments (Karl Pearson, 1894) 

• Method of Maximum Likelihood (R.A. Fisher, 1922) 

• Method of Least Squares (C.F. Gauss, 1810). 

(a) Method of Moments. The k th sample moment is defined by: 



, where n=sample size . 


(245) 


The k tfl population moment is defined by: 

p.' k =\x k f(x)dx . (246) 

The method of moments consists of equating the sample moment and the population moment, 
which means we set m' k -p' k . In general, this leads to k simultaneous (nonlinear) algebraic equations in 
the k unknown population parameters. 

EXAMPLE: Normal distribution: 


and h' 2 =o 2 +H 2 . 

Therefore: 

I x ; .22 

a - — L and — — . 

^ n n n 

From this follows: 


£=*=l5>. 

o 1 = s 2 = -nx 2 )=i £( X t -x ) 2 


(247) 


(248) 


(249) 
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(b) Method of Maximum Likelihood. The likelihood function of a random sample is defined by: 

L(d)=f(x\,X2 ... x n ;0) ■ (250) 

The method of maximum likelihood maximizes the likelihood function with respect to the 

parameter 6. Statistically speaking, this means that the maximum likelihood estimator maximizes the 
probability of obtaining the observed data. 

EXAMPLE: Normal distribution: 


L(n,a 2 )=YlN(x i ;/i,(j 2 ) = 
i=l 


1 


\ n 


a^2n ) 


2a 2 k <X ‘ M)2 


(251) 


Partial differentiation of the “log-likelihood” function £=£n L with respect to jj. and cr 2 yields: 




de = 
do 2 


- rL ~+-l- r (x i -u) 2 = 0 
2cr 2 2cr 4 ' ^ 


Solving for the mean fi and the variance cr 2 we get: 

/i = — 2X- =x and a 2 =—£(x i -x) 2 =s 2 . 
71 * n 1 


(252) 

(253) 


(254) 


It is interesting to observe that the method of moments and the method of maximum likelihood give 
identical estimators for a normally distributed random variable. This is, of course, not so in general. 

Note that: 


• The maximum likelihood estimators are often identical with those obtained by the method of 
moments. 

• The set of simultaneous equations obtained by the method of likelihood is often more difficult to 
solve than the one for the method of moments. 

• Maximum likelihood estimators, in general, have better statistical properties than the estimators 
obtained by the method of moments. 

Invariance Property 


If © is a maximum likelihood estimator of the parameter 0, then g(0) is also a maximum 
likelihood estimator of g(0). For example, if a 2 is a maximum likelihood estimator of the variance a 2 , 
then its square root <T is a maximum likelihood estimator of the standard deviation a\ i.e.. 


< 7 = 




(255) 
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(c) Method of Least Squares. The method of least squares is used in multivariate analyses like 
regression analysis, analysis of variance (ANOVA), response surface methodology, etc. The least-square 
method can be used to estimate the parameters of a two-parameter distribution by transforming the 
cumulative distribution to a straight line, using a probability scale and performing a least-square-curve fit 
to this line. In general, this is a rather unwieldly process unless the cumulative distribution can be obtained 
in closed form such as the Weibull distribution. 

3. Robust Estimation 

Over the last 2 decades statisticians have developed so-called robust estimation techniques. The 
term “robust” was coined by the statistician G.E.P. Box in 1963. In general, it refers to estimators that are 
relatively insensitive to departures from idealized assumption, such as normality, or in case of a small 
number of data points that are far away from the bulk of the data. Such data points are also referred to as 
“outliers.” Original work in this area dealt with the location of a distribution. 

An example of such a robust estimator is the so-called trimmed mean. In this case, the trimming 
portions are 10 or 20 percent from each end of the sample. By trimming of this kind, one removes the 
influence of extreme observations so that the estimator becomes more robust to outliers. (The reader might 
be familiar with the practice of disregarding the highest and the lowest data points for computing athletic 
scores.) A less severe method consists of allocating smaller weights to extreme observations in order to 
reduce the influence of any outliers. This leads to what is called a “Winsorized” mean. Further information 
on this subject matter can be found in Data Analysis and Regression , by Mosteller et al., and in Robust 
Statistics: The Approach Based on Influence Functions by Hampel et al. 

Outliers. At some time or another, every engineer engaged in statistical data analysis will be 
confronted with the problem of outliers. Upon examination of the data, it appears that one or more 
observations are too extreme (high, low, or both) to be consistent with the assumption that all data come 
from the same population. There is no problem when it is possible to trace the anomalous behavior to 
some assignable cause, such as an experimental or clerical error. Unfortunately this is the exception rather 
than the rule. In the majority of the cases one is forced to make a judgment about the outliers, whether or 
not to include them, or whether to make some allowance on some compromise basis. Clearly it is wrong to 
reject an outlier simply because it seems highly unlikely or to conclude that some error in reading an 
instrument must have happened. In fact, in the history of science, attention given to anomalous 
observations has often led to startling new discoveries. For example, the discovery of nuclear fission 
involved enormous energy levels which were originally considered to be outliers. This temptation must be 
strongly resisted. Once one starts handpicking the data, the sample becomes a biased one. The best rule, 
therefore, is to never reject an outlier unless an assignable cause has been positively identified. 

What actually is needed is some objective statistical method to decide whether or not an outlier does 
in fact exist. However, since the outlier problem is by its very nature concerned with the tails of a 
distribution, any statistical method will be very sensitive to departures from the assumed distribution. Due 
to the mathematical difficulties involved, the majority of statistical research has so far dealt with the normal 
distribution. 

One of the earliest methods for dealing with outliers is known as Chauvenet’s criterion (1863). 
This rather strange proposal is stated as follows: If in a sample of size n the probability of the deviation of 
the outlier from the mean is smaller than 1/2 n, the outlier is to be rejected. 

Chauvenet’s criterion is based on the concept of the characteristic extreme (largest or smallest 
value), which was introduced by Richard V. Mises. The characteristic largest value u n is defined as: 

nZ?(u n )=n[l-F(u w )]=l (256) 
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where «=sample size and 


Pr(u<u n )=F(u n ) is the cumulative distribution function. 

Pr(u>u n )=R(u n )=l - F(u n ) . • (257) 

Accordingly, we expect that on the average, one observation in a sample of size n will exceed the value u n . 

Similarly, the characteristic smallest value u\ is defined as: 

nF(«i)=l . (258) 

The number of observations n corresponding to the characteristic largest or smallest value is also called the 
return period T(u n ) or T(u\ ), respectively. This means we expect the maximum or minimum value to 
occur once in a return period. The concept of the return period is sometimes used in engineering design. 

Note that if an event has a probability of p, then the mathematical expectation for this event to 
happen is defined by E-n p. Therefore, if we set E=l, we can make the statement that it takes Up trials on 
the average for the event to happen once. 

It is important to notice that the characteristic extreme value is a population parameter and not a 
random variable. The characteristic largest value is found by solving the above equation for F(u n ): 

F(u n )=l-l/n . (259) 

According to Chauvenet’s criterion, data points will be rejected as outliers if they fall above the critical 
value defined by: 




1 

2 • 


(260) 


In other words, we expect that in a sample of size n, on the average “one-half’ observation will 

exceed the critical value u n c . In the case where data points are rejected, the parameters of the distribution 
are recalculated without the rejected values. If we know the distribution of the population from which the 
sample is drawn, we can determine this critical value from the relation 


F«>=i-y 2n 


(261) 


For a normal distribution N(0, 1) we obtain: 
«50 = 1 “ KoOO 


n = 50 

0 

1 

II 

o 

£ 

“50 — ^0.99 — 2.33 

(262) 

71=500 

O 

O 

II 

h-* 

1 

i 

o 

o 

“500 = 1 - KoOO • 

(263) 


It is seen that when applied to a normal distribution, Chauvenet’s criterion is absolutely 
meaningless. At most, one could use this criterion to flag suspicious data points in a sample. 
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Return periods are often chosen to be equal to twice the service life of a structure when maximum 
loads or winds are fitted to an extreme value distribution (Gumbel distribution). It is seen that the basis for 
the factor two is actually Chauvenet’s criterion. 


C. Sampling Distributions 

The probability distribution of a statistic S is called its sampling distribution. If, for example, the 
particular statistic is used to estimate the mean of a distribution, then it is called the sampling distribution of 
the mean. Similarly we have sampling distributions of the variance, the standard deviation, the median, the 
range, etc. Naturally, each sampling distribution has its own population parameters such as the mean, 
standard deviation, skewness, kurtosis, etc. The standard deviation of a sampling distribution of the 

statistic S is often called its standard error o s . 

1. Sample Mean 

If the random variables X\ , Xi - . .X n constitute a random sample of size n, then the sample mean is 
defined as: 


X = 


IX; 


(264) 


Note that it is common practice to also apply the term statistic to values of the corresponding random 
variable. For instance, to calculate the mean of a set of observed data, we substitute into the formula: 


x = 



(265) 


where the jc, are the observed values of the corresponding random variables. 

First, let us state some general results that hold for arbitrary distributions with mean jl and 
variance o 2 : 


• Infinite population: 


E(X) = p x =p and Var(X)=o% = 



(266) 


• Finite population of size N: 


2 

E(X)=p- x =p and Var(X)=ol=^x^ (267) 

Note that the information embodied in the sample is minimally affected by the proportion that the sample 
bears relative to its population and essentially determined by the sample size itself. In fact, the finite 
population correction factor is usually ignored if the sample size does not exceed 5 percent of the 
population. 
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Large sample size (/V > 30): As a consequence of the Central Limit Theorem, the sampling 
distribution of the mean approaches that of the normal distribution with mean n and variance 

cP-ln. The sampling distribution of the mean (fig. 32) is said to be asymptotically normal. 



FIGURE 32. — Population and sampling distribution. 


• Normal population (a known): If the sample mean is taken from a normal distribution, then 
its sampling distribution is also normal regardless of the sample size. This is due to the 
reproductive property of the normal distribution. 

• Normal population (<r unknown): If the variance of the normal distribution is not known, it 
is required to find a suitable statistic to estimate it. A commonly used estimate is the 
(unbiased) sample variance: 


S 2 YWi-Xp 

n - 1 


(268) 


It can be shown that the statistic: 


S/jn 


(269) 


has die Student-/, or /-distribution with (n- 1) degrees of freedom. This distribution was first obtained by 
William S. Gosset, a 32-year-old research chemist employed by the famous Irish brewer Arthur Guinness, 
who published it in 1908 under the pseudonym “Student.” The /-statistic is also called a /-score. 

Note that because of its smaller MSE, some statisticians prefer to use the (biased) sample variance: 


r2 _l(X i -X) 2 
n 


(270) 


In this case the /-statistic becomes: 


X-// 

S B /in-l 


(271) 
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The Student-? distribution is given by: 


m= 





V 


J 


v+1 

T 


- oo < t < +oo 


where the parameter v=n - 1 is called the “degrees of freedom.” 


(272) 


The variance of the /-distribution for v >2 is v/(v-l) and the kurtosis for v >4 is 3+6/ (v-4). It is 
seen that for large values of v, the variance becomes 1 and the kurtosis becomes 3. In fact it can be shown 
that for v -*», the /-distribution is asymptotically normal (see fig. 33). 



Figure 33. — Student versus normal distribution. 


The /-distribution can be integrated in closed form, but not very easily. Tables are usually given for 
the inverse cumulative distribution; i.e., they present /-values for various percentile levels. 

NOTE OF CAUTION: There exist two categories of /-distribution tables. One type shows percentile 
values tp such that the area under the distribution to the left is equal to p; the other shows values of t a such 

that the area under the distribution to the right; i.e., the tail area, is equal to a. A person who is somewhat 
familiar with the /-distribution should have no difficulty identifying which table to use. 


2. Sample Variance 

Theorem 1. If s 2 is the (unbiased) sample variance of a random sample of size n taken from a 
normal distribution with variance <7 2 , then the statistic: 

x 2 = (n-l)S 2 (273) 

a 1 
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has a distribution with v=n - 1 degrees of freedom whose probability density is: 



0 < X 2 < oo 


(274) 


with mean fl=v and variance <t 2 =2v. 

Note that the % 2 distribution (shown in fig. 34) is a special case of the gamma distribution with the 
parameters «=v/2 and /3=2. 



The distribution of the sample variance s 2 itself is given by: 


f(s 2 ) = 



(275) 


The cumulative % 2 distribution is: 


f (x 2 )= \f(x)dx . 
0 


(276) 


2 2 

Theorem 2. If .Vj and sj are two sample variances of size n i and n 2 , respectively, taken from two 
normal distributions having the same variance, then the statistic: 


F= 


Si 


(277) 
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has a Fisher’s /•’-distribution with: 


vi=n]-l and V2=«2~l degrees of freedom. 


(278) 


The F-distribution is also known as the variance-ratio distribution: 


g(F)=- 


v,+ v 2ir h.fiA 


<V V -2) 
F 2 


Mi Fhr b v 'i 


(v,+v 2 j 

f} 2 


0 <F<« 


(279) 


V 9 

mean H = v _2 f° rv 2 


(280) 


vanance o^- 


2v 2 (V\ +v 2 -2) 

" ^^-2)2^2 -4) 


forv 2 >4 


(281) 


REMARK: The F-distribution is related to the beta function as follows: If X has a beta distribution with 

Vi v 2 . . 

parameters ct= and /3= -y , then the statistic: 


v^l-X) 


(282) 


has an F-distribution with Vi and V 2 degrees of freedom. 


3. Sample Standard Deviation 

For a sample size n>30 taken from an arbitrary population, the sampling distribution of the 
standard deviation S can be approximated by a normal distribution with: 

fi s =a and °* = 2 fc=T) . (283) 

The standard normal random variable is then: 


<x/J2(n-l) 


(284) 
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D. Interval Estimation 


Interval estimates are used to assess the accuracy or precision of estimating a population 
parameter 6. They are defined by the confidence (fiducial) intervals within which we expect a population 
parameter 0 to be located. 

• Two-sided confidence interval: 

P(e L <e<&u)=l-a . (285) 

• One-sided confidence interval: 

- Lower confidence interval: 


P(e L <0)=\-a (286) 

- Upper confidence limit: 

P{e<©u)=\-a . (287) 


The confidence interval is a random variable because the confidence limits (upper and lower endpoints) are 
random variables. 

The probability (1-a ) that the confidence interval contains the true population parameter 6 
is called the degree of confidence or the confidence (level). The corresponding interval is called a 
(1-a) 100-percent confidence interval. The quantity a is called the significance coefficient. 

An analogy given by N.L. Johnson and F.C. Leone in Statistics and Experimental Design states 
that: “A confidence interval and statements concerning it are somewhat like the game of horseshoe tossing. 
The stake is the parameter in question. (It never moves regardless of some sportsmen’s misconceptions.) 
The horseshoe is the confidence interval. If out of 100 tosses of the horseshoe one rings the stake 90 times 
on the average, he has 90-percent assurance (confidence) of ringing the stake. The parameter, just Uke the 
stake, is the constant. At any one toss (or interval estimation) the stake (or parameter) is either enclosed or 
not. We make a probability statement about the variable quantities represented by the positions of the arms 
of the horseshoe.” 

1. Mean 

To obtain a two-sided confidence interval for the mean, we use the /-statistic, which is given by: 


T= 


x-p 

S/4n 


with v=n - 1 


(288) 


and establish the following double inequality 

P\.~ t a/2 < ^ < *a/2 ] = ^ - a 


(289) 


or 
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~ t a/2 < 


X-A* 

S/Jn 


<t. 


a/2 


=l-a 


(290) 


Now we solve the inequality for the desired mean fi and obtain: 


X - t a/2^<f i<X+t a/2-7Z 

V n \in _ 


=l-a with v=n-l 


(291) 


This leads to the following (1-a) 100-percent confidence interval for the mean fl 
(shown in fig. 35): 


s s 

x - t ai 2 -r- <H<x + t a/2 — ~ with v = n-l . 

\in Vn 


(292) 



Following a similar approach as above, we can also establish one-sided upper or lower confidence 
intervals for the mean. For example, the one-sided upper (1-a) 100-percent confidence limit for the mean 
is: 


/j.<x+t a -4= with v=n - 1 . (293) 

a In 

For large sample size (n>30) we simply replace the /-statistic by the z-statistic. Because of the 
central limit theorem, the confidence interval for the mean can also be used for non-normal distributions 
provided n is sufficiently large. 

The most widely used confidence levels 1-a are the 95- and 99-percent levels. For large sample 
sizes, the corresponding two-sided standard scores ( K factors) are zo.025 = l-96 and zo.005 = 2-576. 

If no confidence level is specified, the “one-sigma” interval ( K factor is 1 ) is generally assumed. 
For a two-sided confidence interval, this represents a 68.27 -percent confidence level. 


Error of Estimate. The error of estimate E is defined as the absolute value of the difference 


between the estimator © and the population parameter 6, i.e., E= 


0-0 


. The maximum error of estimate 


is equal to the half-width of the two-sided confidence interval. The maximum error of the mean is, 


therefore: 
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(294) 


E ~ t a/2 h ' 

\n 

Note that the half-width of the two-sided, 50-percent confidence interval for the mean is known as the 
probable error. For large sample size n, the corresponding standard score (K factor) is zo. 25 - 0 - 6745 . 

The above formula can also be used to determine the sample size that is needed to obtain a desired 
accuracy of an estimate of the mean. We solve for the sample size n and get: 


n~ 


t — 

l a/2 £ 


(295) 


2. Variance 

Using the % 2 statistic we can establish a confidence interval for the variance of a normal 
distribution by duplicating the steps we used for the mean (see fig. 36). 



FIGURE 36. — Confidence interval for variance. 


We obtain: 


D (n-l)S 2 2 (n-l)S 2 , 

P<- — ^ — <(T 2 < i -^ r 4 — (=1 -a . (296) 

l %a/2 Xt-a/2 J 

This is an “equal tails” confidence interval. Actually it does not give the shortest confidence interval for d 2 
because of the asymmetry of the x 2 distribution. 

In general, it is desirable that the width of the confidence interval be as small as possible. 
However, the decrease in the interval, which is achieved by searching for the smallest interval, is usually 
not worth the labor involved. 

The confidence of the standard deviation is obtained by simply taking the square root of the above 
formula. It is, of course, also possible to establish upper and lower confidence intervals just as in the case 
of the mean. 
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3. Standard Deviation 


If the sample size is large (n>30), the sampling distribution of the standard deviation S of an 
arbitrary distribution is nearly normal, and we can use its sampling distribution to establish a (1-a) 100- 
percent confidence interval as: 

<cr< s - . (297) 

l+Z a/2 mn~l) l-z a/2 /y2(n-l) 

4. Proportion 

(a) Approximate Confidence Limits. When the sample size n is large, we can construct an 
approximate confidence interval for the binomial parameter p by using the normal approximation to the 

A V 

binomial distribution. To this end we observe that the estimate P=— of the population parameter p has an 
asymptotically normal distribution with mean E(P)=p and variance Var (P)-^. The equivalent 
z-score is: 


P-E( P) _ P~p 

A jVar(F) PI 
n 


where q=l-p ■ 


The two-sided large sample (1-a) 100-percent confidence interval is then given by: 


(298) 


-Z a /2 < 


P-P 

n 


< z a/ 2 • 


(299) 


This is a quadratic equation in the parameter p and could be solved to obtain the approximate 
confidence interval for p. But for the sake of simplicity, we make the further approximation of substituting 
p in the expression for the variance of p appearing in the denominator of the above equation. The 
approximate confidence interval is then determined by: 


P~ z a/2 < P < f >+Z oc/. 2 ^ 


PQ 


=l-a 


(300) 


It is again of interest to find the sample size n that is necessary to attain a desired degree of 
accuracy for the estimate of the binomial parameter p. It is: 


n=p(l-p)\ 



(301) 


Note that if no information is available for the value of p, we set p= 1/2 to obtain the maximum sample 
size. 
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EXAMPLE: Suppose we want to establish a 95-percent confidence that the error does not exceed 
3 percent. Then we must have a sample size of: 



(302) 


(b) Exact Confidence Limits. In reliability work the concern is usually with one-sided 
confidence limits. For a large sample size n one can, of course, establish one-sided limits by the same 
method that was used above for the one-sided confidence limits of the mean. However, it is not very 
difficult to find the exact limits by resorting to available numerical algorithms. These are based upon the 
following mathematical relationships. 

(1) The Upper Confidence Limit. Py is determined from the solution to the following 

equation: 




I 

*=0 


Pu{ l ~PuT k=a ■ 




(303) 


This equation can be solved for the upper limit py by using the equivalency of the binomial sum 
and the beta integral and then equating the beta integral to the /'’-distribution. The result is: 


Pu = 


(* + l )_la (VfiVfe) 


with Vj = 2(x + 1) and v 2 - 2 (n - x) . 


(304) 


equation: 


(n-x) + (x + l) F a (Vj, v 2 ) 

(2) The Lower Confidence Limit. Pi is determined from the solution to the following 


I 

k=x 


n )„k 


P K L( l -PL) n k =a 


(305) 


and solved for the unknown lower limit py in terms of the F-distribution: 


= x+ („- x+ VF a <v vV2 ) with vi=2(n-*+l) and V2=2* . (306) 

NOTE: Comparing these confidence limits for the binomial parameter p with the analogous Bayesian 
confidence limits, it is noticed that the upper limit corresponds to a prior beta distribution with oc= 1 and 
)S=0, whereas the lower limit is obtained by setting these prior parameters to ee=0 and (3=1. 

(3) Two-Sided Confidence Limits. The two-sided confidence limits are obtained by simply 
replacing the significance coefficient a by all for the corresponding upper and lower limits. 


E. Tolerance Limits 

Tolerance limits establish confidence intervals for a desired percentile point where the subscript 
p indicates the proportion to the left of the percentile point of a distribution. They are used in design 
specifications (“design allowables”) and deal with the prediction of a future single observation. The 
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following discussion is restricted to a normal distribution with sample mean x and sample standard 
deviation S. 

1. One-Sided Tolerance Limits 

Here we determine a confidence interval which guarantees a confidence of (1 - a) that at least 
100 p percent of the population lies below (or above) a certain limit. In contrast to the percentile point % p , 
which is a population parameter, the tolerance limit itself is a random variable. 

The one-sided upper tolerance limit (see fig. 37) is defined as: 

F (J =X+K l S . (307) 



FIGURE 37. — One-sided upper tolerance limit. 

Symbolically the above tolerance limit definition can be written as: 

(308) 


where O is the cumulative normal distribution with mean p and standard deviation a. This can be 
expressed somewhat simpler as: 


Pi[F v Z$ p ]= 1 -cc . (309) 

Being a random variable, the tolerance limit has its own probability distribution f(Fy )■ For a 
normal distribution it can be shown that it is a noncentral /-distribution. 

Theoretical work for other distributions is not well developed because of the associated 
mathematical and numerical difficulties. Some approximate tolerance limits have been determined for the 
Weibull distribution. References are listed in the Military Standardization Handbook MIL-HDBK— 5F, 
Vol. 2, 1 November 1990. Another reference on tolerance limits is: Statistical Tolerance Regions: 
Classical and Bayesian, by I. Guttman. 

The one-sided lower tolerance limit is defined as: 

F l = X-K x S (310) 
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with an analogous interpretation for the left-hand tail of the normal distribution. 

Symbolically, we can write: 

Pr[F L <^-p]= l-a . (311) 

Notice that for the normal distribution, we have the relation |j-p = -% p . 

2. Two-Sided Tolerance Limits 


Here we select an interval such that the probability is (l-a) that at least 100 percent of the 
population is contained between the upper limit *U = X+K 2 S and the lower limit F L = X-K 2 S . 

Symbolically expressed, we have: 


Pr 


<P(F u )-^>(F L )>p]=l-a 


(312) 


The mathematical treatment of the two-sided case is surprisingly difficult. The distribution of the 
two-sided tolerance limits (see fig. 38) cannot be expressed in closed form and the numerical effort to 
calculate the exact tolerance limit is quite substantial. For this reason several approximations have been 
worked out and are given in tables. For example, see table 14 in Probability and Statistics for Engineers 
by Miller, Freund, and Johnson, listed in the bibliography. 



NOTE: Two particular tolerance limits have been widely used. Both are using a 95-percent confidence 
level. One is the so-called A-Basic (A-Allowables), which contains 99-percent of the population, and the 
other is the B-Basis (5-Allowables), which contains 90 percent of it. 


F. Hypothesis/Significance Testing 

Here we discuss mainly the Neyman-Pearson theory (1933), also called the classical theory of 
hypothesis testing. The main distinction between hypothesis testing and estimation theory is that in the 
latter we choose a value of a population parameter from a possible set of alternatives. In hypothesis 
testing, a predetermined value or set of values of a population parameter is either accepted or rejected. 
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DEFINITION: A statistical hypothesis is an assumption about a distribution. If the hypothesis completely 
specifies the distribution, it is called a simple hypothesis; otherwise, it is called a composite hypothesis. 

EXAMPLE: Simple H\p=\0 (313) 

Composite H:fJ> 10. (314) 


1. Elements of Hypothesis Testing 

There are four elements of hypothesis testing: 

• Null hypothesis Hq : Baseline or “status quo” claim (design specifications) 

• Alternative hypothesis H \ : Research claim, “burden of proof’ claim 

• Test statistic 0 

• Rejection (critical) region. 

In subjecting the null hypothesis to a test of whether to accept or reject it, we are liable to make two 
types of error: 

Type I error: Reject the null hypothesis Ho when it is true. 

Type II error: Accept the null hypothesis Ho when it is false. 

In the context of quality control of manufactured goods from a certain vendor, the Type I error is 
called the producer’s risk, because it represents the probability of a good product; i.e., one which 
performs according to specifications, being rejected. On the other hand, the Type II error is called the 
consumer’s risk, because it is the probability of a bad product being accepted. 

Symbolically, we have the following conditional probabilities: 

(/?= reject Hq, A=accept Hq) 


oc=Pr (R I H 0 ) and 1 -a=Pr ( A I H 0 ) (315) 

{S=Pr (A I Hi) and l-^=Pr{R\H\) . (316) 

The first column represents wrong decisions and the second column, right decisions. The 
probabilities associated with the right decisions are called the confidence (1— Of) of the test and the power 
(1-/3) of the test. 


In summary. 


• Type I error: 

• Type II error: 

• Confidence: 

• Power: 


Reject good product (a) 
Accept bad product (/3) 
Accept good product (1 - a) 
Reject bad product (1-/3) . 
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NOTE: The probability of a Type I error is also called the significance level. The following two levels are 
most common: 


a= 0 . 05 , called “probably significant” 

0=0.01, called “highly significant.” 

We will, subsequently, notice that the two types of errors have opposite trends; i.e., if one is made 
smaller, the other necessarily increases. Both errors can only be made smaller simultaneously by 
increasing the sample size, thereby increasing the amount of information available. 

To illustrate the general concepts involved in hypothesis testing, we take the example of selecting a 
job applicant for a typing position. As our test statistic 0 , we take the average number x of typing errors 
per page. We define the null hypothesis Hq by the mean number fio of typing errors per page, which 
represents our requirement for the job. The alternative hypothesis H\ is given by the number ji\ of errors 
per page, which we consider unacceptable for the job. 

The next and most important step is to select a criterion based on an actual typing test of the 
applicant that allows us to make a decision whether to accept or reject the applicant. We must categorically 
insist to select this criterion before the experiment is performed. The basic desire is to find an economical 
balance between the two types of errors. This is often not a simple matter because (for a given sample size) 
an attempt to decrease one type of error is accompanied by an increase in the other type of error. In practice 
the “losses” attached to each of the two types of errors have to be carefully considered. One type of error 
may also be more serious than the other. 

In our example we will select a critical value fi c somewhere between (1q and fl\ to define the 
dividing line between the acceptance region and the rejection region. Accordingly, the applicant will be 
accepted if in an actual typing test he or she will have an average typing error jc which is below the critical 
value fl c . Figure 39 illustrates the given example. 



FIGURE 39. — Hypothesis test (Hq: }1=Bq^ H\: 


The different type of errors a and (3 and the associated consequences of making a wrong decision 
should be considered. For example, consider the difference between a civil service job and one in private 
industry in consideration of the prevailing policies pertaining the hiring and firing of an employee and the 
availability of applicants for the job to be filled. 
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NOTE: It is common practice to choose the null hypothesis so that the Type I error becomes more serious 
than the Type II error; in other words, we put the burden of proof on the alternative hypothesis H\ . 

Which value to choose for the significance level depends on the risks and consequences associated with 
committing a Type I error. 

In the following, we present a numerical example of hypothesis testing. 

PROBLEM: A coin is tossed five times. Consider the hypotheses Hq: po=0.5 and H\ : p i=0.7 in which 
p is the probability of obtaining a head. The null hypothesis Hq is rejected if five heads come up. 
Determine the Type I and Type II errors. 

SOLUTION: 


H 0 :po=0.5 H\: p\=0.1 (317) 

a = B(5 1 5,0.5) = Po = 0.5 5 = 0.031=3.1 percent, (318) 

where B represents the Binomial distribution. 

Under the alternative hypothesis, 

4 s 

P= IB(xl5, 0.7) = 1 — pi =1-0.168=0.832=83.2 percent. (319) 

x-0 


Textbooks often cite the analogy that exists between hypothesis testing and the American judicial 
system. The null hypothesis in a criminal court is: 

Hq: “The accused is innocent.” 

The system must decide whether the null hypothesis is to be rejected (finding of “conviction”) or 
accepted (finding of “acquittal”). A Type I error would be to convict an innocent defendant, while a Type 
II error would be to acquit a guilty defendant. Clearly, the system could avoid Type I errors by acquitting 
all defendants. Also Type II errors could be avoided by convicting everyone who is accused. We must 
choose a strategy that steers a course between these two extreme cases. In quantitative terms we must ask 
what kind of “losses” are attached to each kind of error. 

Our judicial system considers a Type I error much more serious: “It is better that 99 percent guilty 
persons should go free than that one innocent person should be convicted” (Benjamin Franklin, 1785). 
The very premise of the system (the accused is innocent until proven guilty beyond a reasonable doubt) 
reflects this concern of avoiding Type I errors, even at the expense of producing Type II errors. 

In the American court of law, the judge often advises the jury to consider the following levels of 
confidence to arrive at a verdict: 

• Preponderance of evidence (> 50 percent) 

• Probable reason (90 percent) 

• Clear and convincing evidence (95 percent) 

• Beyond reasonable doubt (99 or 99.9 percent). 
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NOTE: Many problems people have stem from making Type II errors: 

“It ain’t so much the things we don’t know that get us in trouble. It’s the things we know that ain’t 
so.” Will Rogers (1879-1935). 

“It is better to be ignorant, than to know what ain’t so.” 


2. Operating Characteristic (OC) Curve 

The operating characteristic function (see fig. 40) is the probability of accepting a bad product 
(Type II error) as a function of the population parameter 0: 

L(0)=p L=“ Loss.” (320) 

Notice that when the parameter p=fM) the null hypothesis and the alternative hypothesis coincide. In 
this case, the probability of the confidence and the consumer’s risk are the same, which is to say /3=1 -a. 



FIGURE 40. — Operating characteristic curve. 


Sometimes it is more convenient to work with the power of the test. Remember that the power is 
the probability of rejecting the null hypothesis when the alternative is true. The so-called power function 

7t(d) of the test is related to the operating characteristic function by: 

tt(0)=1-L(0) . (321) 

Any test becomes more powerful as the sample size n increases. 

3. Significance Testing 

It is often difficult to specify a particularly meaningful alternative to the null hypothesis. It might 
also not be feasible in terms of resources and economy to determine the Type II error. In this case the 
alternative hypothesis Hi is composite ( i.e., ). These tests are also called null hypothesis 

tests. 
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If the test statistic© falls in the acceptance region, we are then reluctant to accept the null 
hypothesis because the protection against a Type II error (consumer’s risk) is not known. Therefore, one 
usually prefers to say that the null hypothesis “cannot be rejected” or we have “failed to reject” the null 
hypothesis. 

There are two different classes of hypothesis tests, depending on the type of rejection criterion that 
is used. 

(a) One-Sided Criterion (Test). The null hypothesis is rejected if the value of the test statistic 0 
lies in one (upper or lower) tail of the sampling distribution associated with the null hypothesis Hq (see 
fig. 41). This is also called a one-tailed test. 

In case of a significance test, the alternative hypothesis H\ is called one-sided : 

H Q :0=e 0 ,H l :d>e 0 ,orH l :6<d 0 . (322) 


H 0 H, 



(b) Two-Sided Criterion (Test). Here the null hypothesis is rejected if the value of the test 
statistic © falls into either tail of the sampling distribution of the associated null hypothesis Ho (see fig. 
42). 


H 0 Hi 



FIGURE 42. — Two-sided hypothesis test. 

In case of a significance test, the alternative hypothesis Hj is called two-sided : 

Hq.6=0q H v G*0 0 . (323) 
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Figure 43 illustrates the decision steps taken in a significance test. 



(Cannot Reject H 0 ) 

FIGURE 43. — Significance test (a=0.05) 
EXAMPLE: Significance testing concerning one mean: 


Test statistic: 


x-ll n 

t = — £2- , v=n-l 
s An 


Ho : IU=W) 


(324) 

(325) 


H o 

Reject Hq 

0 

A 

1 

ft 

o 

t > t a 

o 

t > taj 2 or t < -tall 


REMARK: Tests concerning other statistics can be found in most statistics reference books. 


G. Curve Fitting, Regression, and Correlation 


1. Regression Analysis 

Regression assumes the existence of a functional relationship between one dependent random 
variable y and one or more independent nonrandom variables, as well as several unknown parameters: 

Y = f(x\, X 2 . . ■ x n ; 0i, 62 . . . d m ) + E (326) 

Y = Response (criterion) variable 

Xi = Prediction (controllable) variable 

Oj = Regression parameters 
E = Random error. 
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The regression analysis consists of estimating the unknown regression parameters Oj for the 
purpose of predicting the response Y. Estimation of the (population) regression parameters is done by the 
method of least squares, which is almost universally chosen for “curve fitting.” Accordingly, we 
minimize: 

S = ~f(*V x 2 • • x n ’ e v 6 2 - e m (327) 

by setting the partial derivatives with respect to the m unknown parameters Oj equal to zero. The resulting 
m simultaneous equations are called normal equations. They are, in general, nonlinear and difficult to 
solve unless the parameters enter linearly. Sometimes a suitable transformation can be found which 
“linearizes” the problem. Another approach of linearization is to find good initial estimates for the 
parameter, so that the nonlinear functions can be approximated by the linear terms of the series expansion. 

2. Linear Regression 

Linear regression is expressed as: Y=a + fix + E. 

We estimate the regression parameters a and ft by the method of least squares (see fig. 44). 


e, = Vertical deviation 
of a point y, from 



Figure 44. — Linear regression line. 


Setting S='£ef = Y. yj , we solve J^=0 and ^=0. 

(a) The Normal Equations are given as: 


(328) 


an+P? t x i ='Zy i 


( 329 ) 



(b) Abbreviations. We next introduce the following standard notation: 

S xx = l(x i -.*) 2 = Xx?-(lx i ) 2 /n 

s yy = If 3V ~y ) 2 = I yf -(l y t ) 2 /« 

^ = If*, -xXy t ~y) = 1 x t y t - (!>,. )(ly ( > • 

(c) Least-Square Estimators. The least-square estimators are: 

a = y- fix and ft = — — with E(a) = a and £(/?) = /J. 

•A rr 


(330) 


(331) 


3. Gauss Markov Theorem 

The Gauss Markov theorem states that: 

• The estimators cc and /3 are unbiased. 

• The estimators a. and ft have the smallest variance, i.e., they are the most efficient ones. 

This theorem is independent of the probability distribution of the random variable y. The variance 
O' 2 of the random error E is usually estimated by: 

° 2=s l =^2^-(« + ^i)] 2 • (332) 

The division n- 2 is used to make the estimator for a 2 unbiased. The term s e is called the standard 
error of the estimate. 

The computational form of o is: 


gi- ^yy. 

~ n-2 


(333) 


4. Normal Regression Model 

The normal regression model is expressed as: 


/(» I x t ) = — ]= exp - 1 
o^ln 2 


y, -(« + /&, ) 


for -oo <y. <oo . 


(334) 


The maximum likelihood estimators for a, /?, and are identical to the least square estimators 


r\ 

except that o has the divisor n instead of (n-2). 
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The following random variables are used to establish confidence intervals and to perform 
significance tests concerning the estimated regression parameters a and f. 


Intercept or. 


t = 


a-a nS xx 


s e ^S^+x 


v 2 


v=n-2 


(335) 


Slope (3: 


t=&- 

i p 


S e V' 5 xx 


v=n-2 


(336) 


EXAMPLE: Slope /?: 


Ho '- P~Po 


(337) 


Alternative hypothesis H\ 

P<Po 

P>Po 

P±Po 


Reject Hq if 

t < - t a 
t>t a 

t < -tail or * > toj 2 


Note that if the independent variable is time, the regression line is also called a trend line. 
Confidence interval about the regression line: 


y- tail s ej~ + ** - * ■- <y<y+t a /2 Jr + 


jl , {X-X ) 2 


1" S x 


\ n S 


XX 


= 1 -a (v = n-2) . (338) 


Prediction interval about the regression line: 


„ L i (x-jc ) 2 

y-tal 2 s ejl + _ + — 7 , < y < y + t a /2 s e 

V n S xx 


1 i (x-x y 
i+-+- - 


n 


= 1 -a (v = n - 2) .(339) 


The minimum width occurs at x = x (see fig. 45). 
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Note that the width of the prediction interval does not approach zero as n approaches infinity, 
expressing the inherent uncertainty of a future observation. Also, the intervals become increasingly wide 
for x — values outside the range of data. 

5. Nonintercept Linear Mode 

Sometimes both theoretical considerations or empirical data analysis indicate that the linear 

regression line goes through zero; i.e., we impose the condition that a=0. The regression line (see fig. 46) 
is then simply given by: 


yi-bxi . 


(340) 



The regression parameter b is again determined by minimizing the error sum of squares: 

Minimizing S= foc f ) 2 we set ^§-=0 to obtain 

1 1 an V 


I xf 


(341) 


The parameter b is again unbiased, such that E(b)=f3 and its variance is a % = 


2 _ cr 


5 >- 


2 ' 
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The confidence interval for the regression parameter /3 can then be established from the fact that 
the following statistic follows a /-distribution: 


/ = - — — \ I with v=n - 1 degrees of freedom 


s e A " 


and 


s2 e =~ii^(yr bx i) 2 ■ 

The two-sided prediction interval can be shown to be: 


y=bx±t a/2 s e j(l+x 2 /£x?'j with v=n- 1 


(342) 


(343) 


(344) 


The one-sided prediction interval is obtained by replacing tajj by t a and choosing the proper sign 
for the upper or lower limit. 

REMARK: If it is difficult to decide on physical grounds whether the model should contain an intercept 
term, it is common practice to initially fit an intercept model and perform a nullity test for the intercept term 

(i.e., test the null hypothesis Ho-GC=0). If this term is not significant, then the model is refitted without the 
intercept term. 

6. Correlation Analysis 

It is assumed that the data points (jq, y,) are the values of a pair of random variables with joint 
probability density: 


fix, y)=g(y I x ) h\{x) (345) 

where: 

g(y \x)=N(a+ftx,a 2 ) 

tyx^N^o?) (346) 

h 2 (y)=N(ti 2 ,ol) 


f(x,y)= 


1 

2na x a 



(347) 


It can be shown that: 


H 2 = a + At/j and a\ = a 1 + 


(348) 


Define: 


p- 1-^t or p= 1-^j ■ 

\ e\ 


(349) 
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Also: 


P=^~P and o 2 =ol(\-p 2 ) 


(350) 


Then 


exp ~W%\ 

f(x,y)= ^ which is the Bivariate Normal Distribution, 

2na { a 2 ^-p 2 


(351) 


where 


Q(x)= 


\x-m~ 

2 

—2p 

-x-H 

~y~H 

+ 

~y-p 2 

2 ' 

[ a i J 


L a i J 

°2 \ 


L a 2 J 



(352) 


7. Other Regression Models 

Three other regression models are polynomial, multiple, and exponential. 

Polynomial regression: 

y = P 0 +P l x+p 2 x 2 + ...B p xP+£ . (353) 

Multiple regression: 

y = P 0 + Pl x l+p2 x 2 + ‘-- B r x r +e • ( 354 ) 

These two models are still linear regression models because they are linear in the unknown regression 
parameters. 

Using matrix notation, we can write: 

Y=X$ + e (355) 

S = e T e = (Y-XP) T (Y-XP). (356) 

Minimizing S we obtain 

j- = 2X T (Y-XP) = 0 (357) 

and 

P = (x T Xy i X T Y. (358) 

The matrix C=(X7 X) -1 X T is called the pseudo-inverse of X. 
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For example, if we consider the exponential model: 


y = a e 


fix 


We “linearize” it by taking logarithims of both sides, to obtain: 

£ny = £na+Px. 


(359) 


(360) 


8. Sample Correlation Coefficient 

Sample correlation coefficient scattergrams are shown in figure 47. The correlation coefficient 
allows a quantitative measure of how well a curve describes the functional relationship between the 
dependent and independent variables. 



s 2 =X(y,-y /) 2 

Unexplained Variation 


• — ^ y 

: • * Sl=Z(y r y) 2 

^ Total Variation 

x 

FIGURE 47. — Sample correlation coefficient (scattergrams). 

Decomposition of the Total Sum of Squares . The derivation of the correlation coefficient is 
based on the fact that the total sum of squares of the deviation of the individual observations y\(i=\,2 ... 
n) from the mean y can be decomposed in two parts if the regression model is a polynomial of the form: 

y = bQ+b^ x +£>2 x^+...b k x k . (361) 

The error sum of squares is then given as: 


S=le} = hy r y,f 

1=1 1=1 

n 2 

= X [y t - ( b 0 + i b x x +, b 2 X 2 + . . .b k x k )] . 


(362) 
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Partial differentiation with respect to the regression parameters /?, yields the normal equations: 


^-2l[y i -(b 0+ b l x + ...b k x^0^e i =0 

^=-2l[y i -(b 0 +b l x+...b k x k j\x i =0=*'Ze i x i =0 (363) 

= -2l[y ( - -(b 0 +b l x+...b k x k j\x k =0=*Ze i xj‘ =0 . 

Therefore: 

'L£i9i='L£j(b 0 +b ] x+...b k x k )='Ze i +'Ze i x i +...'Z£ i xf = 0. (364) 

This is a significant result because it reveals that the correlation between the error e and the estimate y is 
zero, or one can also say that all the information contained in the data set y has been removed. 

The total sum of squares denoted by SST can then be written as: 

sst=s„ = i(yryf=z(yry + yryi) 2 

= 'L(yi-y) 2 +'L(y i -y i ) 2 +2J d (y i -y)(y i -y i )- (365) 

The last term on the right-hand side can be seen to be zero as follows: 

!(?,- - J^i “>,•)= l(9j ~y)£i = 19 i £i -yl£i = o . (366) 


Therefore, we have the final result that the total variation can be decomposed as follows: 


SST = I(y i -y) 2 =I(y J .-y) 2 +S(y ; -y / ) 2 =SSR+SSE (367) 

where SSR is the regression (“explained”) sum of squares and SSE the error (“unexplained”) sum of 
squares. 

The coefficient of determination is the ratio of the explained variation to the total variation. Since 
it is always positive, we denote it by r 2 : 


2 _ - y ) 2 Xty; -y ) 2 _ s 2 

r ~l(yi-y) 2= "I (*-y) 2= ~ s} * 

Note that r 2 does not depend on the units employed because of its nondimensional nature. 


(368) 
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The positive square root is called the (nonlinear) correlation coefficient: 


r = 


! I(>v-.y ) 2 

1 — o' . 

\ l(yi-y ) 2 


( 369 ) 


The correlation coefficient measures the goodness of fit of the curve. For instance, if y ( - ~ y h then 

r a* 1 . When r = 0, the relationship between the variables is poor with respect to the assumed regression 
model. 

NOTE: The correlation coefficient as defined in equation (3) above can become imaginary for nonlinear 
regression when the assumed model does not fit very well. However, the upper limit is always 1 . 

9. Linear Correlation Coefficient 


Here it is assumed that the relationship between the two variables is linear. Therefore: 


S 2 ='L(y i -y) 2 where y ( - =a+bx i 
S 2 = X>f -2 'L(a+bx i )y i + 'Z(a+bx i ) 2 . 


(370) 


Remember: 


a-y-bx 



(371) 


S 2 = £y 2 -2 a'Zy i -2 b'^x i y i +na 2 +2ab'£x i +b 2 X* 2 
S 2 =X)f -2(y-bx)'£y i -2b'Zx i y i +n(y-bx) 2 +2(y-bx)b'£x i +b 2 '£ t xf 
S 2 = 'Ly 2 -2y Xy,- +2b*Xy,- -2 frX*j y,- +ny 2 -2nbyx+nb 2 x 2 
+2byJ i x i -2b 1 xJ d x i +b 2 X*? 


S 2 =[Lyf-ny 2 ]-2b[Zx i y i -nxy]+b 2 [£x?-nx 2 ] 


Recall S 2 = S XX so that 


(372) 
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(376) 



r = 




\ $xx$yy 

This is called the hnear correlation coefficient. The linear correlation coefficient can also be written as: 

I (Xj-xXyj-y) 


V l(x i -x) 2 'Z(y i -y) 2 

Usually the term correlation coefficient is used to mean linear correlation coefficient. 


(378) 


The sample correlation coefficient r = p is a biased estimator of p in the bivariate normal 
distribution. 


NOTE: The method of maximum likelihood applied to the parameter 0 \ , < 72 , and p of the bivariate normal 
distribution yields the same estimator for the correlation coefficient as above. 

The correlation coefficient can assume positive and negative values (-1 < r < 1), depending on the 
slope of the hnear regression Une (see fig. 48). 




FIGURE 48. — Positive versus negative correlations. 
From equation (368) we can write: 



(379) 


—Sj (l-r 2 ) . (380) 

The variance s 2 is, so to speak, the “conditional” variance of y given x. 

For r=± 1, the conditional variance j 2 =0; i.e., the data points fall on a straight line (perfect 
correlation). 

Sometimes equation (379) is written in the form: 
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x 1 00 percent . 


(381) 


2 

S 2 2 


Thus, 100 r 2 is the percentage of the total variation of the dependent variable of the dependent variable y 
which is attributed to the relationship with x. This is also true for the nonlinear case. 

EXAMPLE: If r= 0.5, then 25 percent of the variation of y is due to the functional relationship with x. Or 
we might say that a correlation of r=0.6 is “nine times as strong” as a correlation of r=0.2. 

NOTE: There are several serious pitfalls in the interpretation of the correlation coefficient. It has often been 
said that it is the most abused statistical quantity. 

PITFALL 1 : The linear correlation coefficient is an estimate of the strength of the linear association 
between the random variables. See figure 49 for a quadratic relationship with zero correlation. 



FIGURE 49. — Quadratic relationship with zero correlation. 

PITFALL 2: A significant correlation does not necessarily imply a causal relationship between the two 
random variables. For example, there may be a high correlation between the number of books sold each 
year and the crime rate in the United States, but crime is not caused by reading books (spurious 
correlation). 

Often two variables are having a mutual relationship with a third variable (e.g., population size) 
which produces the correlation. These cases can sometimes be treated by “partial” correlation. 


10. Sampling Distributions of r 

To establish confidence intervals and to perform significance tests, one would need to know the 
probability distribution of r for random samples from a bivariate normal population. This distribution is 
rather complicated. However, R.A. Fisher (1921) found a remarkable transformation of r which 
approximates the normal distribution. 

11. Fisher Z-Transformation 

The Fisher Z-transformation is given as: 

z = ±£n^=tanh- l r (382) 
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with 




= li n l±P 

Z 2 m \-p 


and <7 


2 _ 1 
2-^=3 


for n > 30 


EXAMPLE: 


r= 0.70 
Z= 0.867 


CONFIDENCE INTERVAL: 


n=30 

a=0.05 


Z— 


"an 

v'rt-3 


<n 7 <z+- 


g/2 

Jn-3 


(383) 


(384) 

(385) 


(386) 


where Zq/ 2 = 1-96 for 95 -percent confidence: 


0. 867-^M < u < 0. 867+ W 
\ 27 ^ V27 

0.490 <p z < 1.244 

or for the correlation coefficient: 

0.45 < p < 0.85 . 

Estimates of correlation coefficients based on small sample sizes are not very reliable. For 
instance, if the sample size n- 20, the calculated correlation coefficient has to exceed the critical value 
r^O.377 to be significant at a=0.05. 


(387) 

(388) 

(389) 


H. Goodness-of-Fit Tests 

Goodness-of-fit tests examine how well a sample of data agrees with a proposed probability 
distribution. Using the language of hypothesis testing, the null hypothesis Hq is that the given data set 
follows a specified distribution F(x). In most applications the alternative hypothesis H\ is very vague and 
simply declares that Hq is wrong. In other words, we are dealing with significance tests in which the Type 
II error or the power of the test is not known. Moreover, in contrast to the usual significance tests where 
we usually look for the rejection of the null hypothesis in order to prove a research claim, goodness-of-fit 
tests actually are performed in the hope of accepting the null hypothesis Hq. 

To make up for the lack of existence of a well-defined alternative hypothesis, many different 
concepts have been proposed and investigated to compare the power of different tests. The result of this 
quite extensive research indicates that uniquely “best” tests do not exist. There are essentially two different 
approaches to the problem of selecting a proper statistical model for a data set — probability plotting 
(graphical analysis) and statistical tests. 
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1. Graphical Analysis 

Graphical techniques can be easily implemented and most statistical software packages furnish 
specific routines for probability plotting. They are very valuable in exploratory data analysis and in 
combination with formal statistical tests. In the former, the objective is to discover particular features of the 
underlying distribution, such as outliers, skewness, or kurtosis (i.e., thickness of the tails of the 
distribution). The old saying that a picture is worth a thousand words is especially appropriate for this type 
of analysis. 

Of particular importance is the so-called empirical distribution function (E.D.F.) or sample 
cumulative distribution. It is a step function that can be generally defined by: 

F ’ <x ‘> = li=E+l foro£cSI < 390 > 

with the observed ordered observations x\ < X 2 ... < x n . Several values for the constant c are in vogue for 
the plotting (rank) position of the E.D.F. The “midpoint” plotting position (c=0.5) has been found to be 
acceptable for a wide variety of distributions and sample sizes. The “mean” plotting position (c=0) is also 
often used. Another one is the so-called “median” rank position which is well approximated by c= 0.3 
(Benard and Bos-Levenbach, 1953). In practice, it does not make much difference which plotting position 
is used, considering the statistical fluctuation of the data. However, one should consistently use one 
particular plotting position when comparing different samples and goodness-of-fit tests. 

A major problem in deciding visually whether an E.D.F. is conforming to the hypothesized 
distribution is caused by the curvature of the ogive. Therefore, it is common practice to “straighten out” the 
plot by transforming the vertical scale of the E.D.F. plot such that it will produce a straight line for the 
distribution under consideration. A probability plot is a plot of: 


Zi=G-HF n (xi )) ( 391 ) 

where G -1 (.) is the inverse cumulative distribution. Sometimes the z-score is placed on the horizontal axis 
and the observations x ; on the vertical axis. By proper labeling of the transformed scale, one obtains the 
corresponding probability graph paper which is commercially available or can be generated with the 
computer graphics package. 


2. x 2 Test 

This test is the oldest and most commonly used procedure for examining whether a set of 

observations follows a specified distribution. The theoretical work underlying the x 2 test was done in 
1875 by the German physicist Friedrich Helmert. The English statistician Karl Pearson (1857-1936) 
demonstrated its application as a goodness-of-fit test in 1900. The major advantage of the test is its 
versatility. It can be applied for both discrete and continuous distributions without having to know the 
population parameters. The major drawback is that the data have to be grouped into a frequency 
distribution (histogram), which requires a fairly large number of observations. It is also usually less 
powerful than tests based on the EDF or other special purpose goodness-of-fit tests. The test statistic uses 
the observed class frequencies O, of the histogram and the expected theoretical class frequencies £,, which 
are calculated from the distribution under consideration, and is defined by: 




E i 


k 0} 

= — n with E- > 5 

i= 1 E i 


(392) 
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where k is the number of class intervals of the histogram and n the total number of observations. The 
following constraint exists between the observed and the expected frequencies: 

£0. =££.=„. (393) 

The sampling distribution of this statistic is approximately the x 2 distribution with degrees of 

freedom v=k- 1 -m where m is the number of population parameters that have to be estimated from the 
data. For this approximation to be valid, the number of expected frequencies £, should be greater than 
five. If this number is smaller, then adjacent class intervals can be combined (“pooled”). Class intervals do 
not have to be of equal size. 

The x 2 test statistic is intuitively appealing in a sense that it becomes larger with increasing 
discrepancy between the observed and the expected frequencies. Obviously, if the discrepancy exceeds 
some critical value we have to reject the null hypothesis Ho which asserts that the data follow a 
hypothesized distribution. Since the alternative hypothesis is not specified, this procedure is also called the 

X 2 test of significance. Typical significance levels are a=0.10, 0.05, and 0.01. When selecting a 
significance level it is important to keep in mind that a low level of significance is associated with a high 
Type II error, i.e., the probability of assuming that the data follow a suggested distribution when they are 
really not. Therefore the higher we choose the significance level, the more severe is the test. This aspect of 
the goodness-of-fit test is, at first sight, counterintuitive, because, as was mentioned above, we are usually 
interested in rejecting the null hypothesis rather than in accepting it. 

Table 5 illustrates the procedure of applying the x 2 test. It shows a table of the observed and 
expected frequencies of tossing a die 120 times. 

TABLE 5 . — Procedure of applying the x 2 test. 

Die face 1 2 3 4 5 6 

Observed frequency 25 17 15 23 24 16 

Expected frequency 20 20 20 20 20 20 

The test statistic yields ^ 2 =5.00. Since the number of class intervals is 6 and no population 

parameter is estimated, the degrees of freedom are v=6 - 1=5. If we choose a significance level of a=0.05 

the critical value is J 2 o.05=H-l- Therefore, we accept the null hypothesis, which means we assume that 
the die is fair. 

We must also look with skepticism upon a x 2 value that is unreasonably close to zero. Because of 
the natural statistical fluctuation of die data, we should not expect the agreement between the observed and 
the expected frequencies to be too good. 

The problem of a small x 2 value is illustrated by the strange case of monk Gregor Mendel’s peas. 
Writing in 1936, the famous English statistician R.A. Fisher wondered if, given the statistical fluctuations 
of experimental data in the field of genetics (Mendel, 1822-1884), Mendel’s results were too good. In 

effect, Fisher tested the left-hand side of the x 2 distribution for a too low x 2 value. Exami nin g Mendel’s 
data conducted over an 8-year interval, he found the x 2 value to be £ 2 =4 1.606 with 84 degrees of 
freedom. This corresponds to a level of significance of a=2.86 x 10— 5 . Therefore one would expect to 
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find such a low x 2 value only three times in 100,000 experiments. Perhaps Mendel was deceived by an 
overly loyal assistant who knew what results were desired. 

It has been said by some critics of the x 2 test, that it will always reject the null hypothesis if the 
sample size is large enough. This must not necessarily be the case. 


3. Kolmogorov-Smirnov Test 

The Kolmogorov-Smirnov test (see fig. 50) is one of many alternatives to the x 2 test. The test is 
called an E.D.F. test because it compares the empirical distribution function with the theoretical cumulative 
distribution. The E.D.F. for this test is defined as: 


F n (xj)=i/n (394) 

where the jc,’s are the ordered n observations. With each increase of an x value, the step function takes a 
step of height 1/n. The test statistic is: 


D - max I F n (jc/) - F{x) I . (395) 

If D exceeds a critical value D a , the null hypothesis is rejected. 

1 

0.9 

F(x) 0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 

0 *1 *2 x n ' * 

FIGURE 50. — Kolmogorov-Smirnov test. 

The Kolmogorov-Smirnov test can also be used for discrete distributions. However, if one uses the tables 
which assume that the distribution is continuous, the test turns out to be conservative in the sense that if 
H 0 is rejected, we can have greater confidence in that decision. In cases where the parameters must be 
estimated, the Kolmogorov-Smirnov test is not applicable. Special tables exist for some particular 
distributions such as the exponential and the normal distribution. For a normal distribution for which the 
mean and variance have to be estimated, the following asymptotic critical D values (table 6) can be used if 
the sample size n > 20. 


a 

D c 


TABLE 6. — Normal distribution. 

0.01 0.05 0.10 

1.031 0.886 0.805 

~Jn 4n 4n 
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Other E.D.F. statistics that measure the difference between F n (jc) and F(x) are found in the literature. A 
widely used one is the quadratic statistic: 


Q= n I { F n( x )~ F ( x ^\ L Y(x)dF(x) (396) 


where \|r(x) is a suitable weighting function for the squared differences under the integral. When y/(x)=l, 

the statistic is the Cramer von Mises statistic, now usually called W 2 , and when i//(;t)=[F(jt){ 1 -F(jc)] -1 

the statistic is the Anderson-Darling (1954) statistic, called A 2 . The E.D.F. test associated with the latter 
statistic is recommended in Volume 2 of MIL-HDBK-5F for testing the normality of data. 


I. Quality Control 


In recent years engineers have witnessed a dramatic revival of interest in the area of statistical 
quality control (SQC). One of its new features is a “process orientation” according to which emphasis is 
shifted towards the improvement of a product during the engineering design and manufacturing phases 
rather than attempting to control quality by inspecting a product after it has been manufactured. Time- 
honored maxims are being quoted and rediscovered. Some of them are: “It is more efficient to do it right 
the first time than to inspect it later.” “One should fix processes and not products.” “All or no inspection is 
optimal.” The oldest and most widely known is, of course, that “one cannot inspect quality into a 
product.’’' 

In addition, the engineering design phase is being modernized and improved by recognizing the 
importance of applying the concepts of experimental design, which had been miserably neglected or 
overlooked in the field of engineering in the past. This is now often referred to as off-line quality control. 
Someone has recently remarked that if engineers would have a better background in engineering 
probability and statistics, it would not have been possible for Professor Genichi Taguchi to dazzle them 
with what are essentially elementary and simple concepts in the design of experiments that had been 
known for a long time by practicing statisticians. And last, but not least, a new element has been 
introduced called “the voice of the customer” (VOC), which is the attempt to get some feedback from the 
consumer about the product. 

SQC centers primarily around three special techniques: (1) determining tolerance limits in the 
design phase, (2) setting up control charts for on-line quality control or statistical process control (SPC), 
and (3) devising acceptance sampling plans after the product has been manufactured. 

1. Acceptance Sampling 

Acceptance sampling is designed to draw a statistical inference about the quality of a lot based on 
the testing of a few randomly selected parts. While acceptance sampling is usually considered part of 
quality control, it is very important to understand that it hardly exercises direct leverage over process 
control. Because many contracts still contain requirements for submitting acceptance sampling plans, it is 
desirable for an engineer to know the basic underlying statistical concepts. 

Statistical methods of acceptance sampling are well developed for a variety of samp lin g procedures 
such as single, multiple, and sequential sampling. A single sampling plan consists of drawing a random 
sample and develops criteria for accepting or rejecting the lot. In a multiple sampling plan, the sampling 
occurs in several stages. At each stage a decision is made whether to accept or reject the lot or whether to 
continue taking another sample. The concept of multiple sampling can be generalized to sequential 
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sampling, in which a decision is made to accept, reject, or continue sampling after each observation. It has 
turned out, against original expectations, that sequential sampling can significantly reduce the amount of 
inspection. 

Acceptance sampling is further classified, depending on whether it applies to attribute scales or 
measurement (variable) scales. In general, inspection by attributes is less expensive than by 
measurements. However, inspection by variables provides more information than attributes and can be 
better used to improve quality control. 

To facilitate the design and use of acceptance sampling plans, several standard plans have been 
published. Among the most widely used are the Military Standard 105D Tables (1963) for attribute 
sampling and the Military Standard 414 Tables (1957) for measurement sampling. 

The subsequent discussion will be limited to single sampling plans. 

• n=Sample size 

• N= Lot size 

• x=Number of defective items in sample 

• c=Acceptance number (accept lot when x < c) 

• Ho : po=Acceptable quality level (AQL) denoting a “good” lot 

• Hi : pi=Lot tolerance percent defective (LTPD) denoting a “bad” lot 

• oc=Pt (Reject Ho I po)=P r °ducer’ s risk=probability of rejecting a good lot 

• /J=Pr (Accept Ho I pi)=consumer’s risk=probability of accepting a bad lot. 

A given sampling plan is best described by its OC curve, which is the probability of the 
consumer’s risk for each lot proportion defective p. This probability can be calculated using the 
hypergeometric distribution as follows: 


h(x)= 



(397) 


where Z)=A),=defective items in the lot. 

The cumulative hypergeometric distribution yields the probability of accepting a lot containing the 
proportion of defectives p: 


c 


L(p)= Y d h(x\n,N,D). 

x=Q 


(398) 


For large lot sizes it is acceptable to approximate the hypergeometric distribution by the binomial 
distribution. A sketch of the OC curve is given in figure 51. The OC curve always passes through the two 

points {po , oc) and (p\, (5), which are reached by agreement between the consumer and producer. 
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FIGURE 5 1 . — OC curve for a single sampling plan. 


Sometimes rejected lots are screened (100-percent inspection) and all defective items are replaced 
by good ones. This procedure is called rectifying inspection. In this case the sampling plan is described 
by its average outgoing quality or AOQ curve. The formula for it is derived under the above assumption 
that all defectives in a rejected lot are replaced by good items before their fined acceptance. Accordingly we 
have: 


AO Q=p L(p) + 0 x (1 -L(p)) 


(399) 


or 


AOQ=p L(p) . (400) 

In general, there will be a maximum average outgoing quality as a function of the incoming lot 
quality p. This maximum is called the average outgoing quality limit (AOQL). Thus, no matter how bad 
the incoming quality becomes, the average outgoing quality will never be worse than the AOQL. 

2. Control Charts 

While studying process data in the 1920’s, Dr. Walter Shewhart of Bell Laboratories first 
developed the concept of a control chart. Control charts can be divided into control charts for 
measurements (variables) and for attributes. Measurements are usually continuous random variables of a 
product characteristic such as temperature, weight, length, width, etc. Attributes are simply judgments as 
to whether a part is good or bad. 

All control charts have the following two primary functions: 

• To determine whether the manufacturing process is operating in a state of statistical 
control in which statistical variations are strictly random. 

• To detect the presence of serious deviations from the intrinsic statistical fluctuations, called 
assignable variables icauses), so that corrective action can be initiated to bring the process 
back in control. 

Control charts are defined by the central line , which designates the expected quality of the 
process, and the upper and lower control limits whose exceedance indicates that the process is out of 
statistical control. However, even when the sample point falls within the control limits a trained operator is 
constantly monitoring the process for unnatural patterns such as trends, cycles, stratification, etc. This 
aspect is called control chart pattern analysis. 
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The following lists the three types of control charts for measurements: 
X -Chart: 

Given: fi, a Central Line: fi 

UCL: ji + A a where A=3/ -Jn 
LCL: n~Ao 

Given: x, s Central Line: x 

UCL: x + A i s 

= j* 

with x=-r'£x i 
k i = 1 

LCL: x -A\ s 

In on-hne quality control, the biased sample standard s is used: 


(401) 


s= 


l'L(x i -x) 2 
V n 


(402) 


Typical sample sizes n are 4, 5, 6, or 7. These relatively small sample sizes often come from the 
selection of subgroups. The control factor A\ is obtained from the expected value of the standard deviation 

s, which is fls= c 2 <7 where: 


c 2 — 




where T (x) is the gamma function. 


(403) 


Therefore, the control factor is A \ =A/c 2 = 3/^c 2 v n ) . 

This and subsequent control factors are calculated under the assumption that the measurements 
come from a normal population. 


Given: x, R 


Central Line: jc 

= - _ i k 

UCL: *+ A 2 r with R=f£K. 

K i = 1 


(404) 


LCL: x-A 2 r 


110 



The range R is widely used in quality control because it can be easily obtained if the sample size is 
small. However, since R depends only on two observations, namely the maximum and the minimum value 
of the sample, it contains less information than the standard deviation. This loss of information is 
acceptable for a small sample size. A practical rule is to use the standard deviation s when the sample size 
is larger than 12. 

In general, the probability distribution of the range R cannot be expressed in simple explicit form. 
For the normal case the first four moments of the distribution of the range R have been extensively 
tabulated as a function of the sample size n. For the control charts we need only the mean and standard 

deviation of the range R. They are denoted by HR=d 2 <7 and ( 7 R=d^ a. The above control factor is then 
given by A 2 =A/d 2 - 

R-Chart: (n < 12) 

Given: a Central Line: flR=d 2 <7 (405) 

UCL: D 2 o 
LCL: D\ a 

where T> 2 = dl + 3^3 and X)\=d 2 - 3^3 . (406) 

Given: R Central Line: R 

UCL: D A R 
LCL: £> 3 R 

where D 4 = 1 ) 2 / d 2 = (d 2 + 2 >d^)/d 2 and Dy=D\ld 2 ={d 2 ~$d'$)ld 2 - (407) 

s-Chart: (« > 12) 

Given: <7 Central Line: Hs=C 2 0(411) (408) 

UCL: B 2 a 


LCL: B\ a . 

The control factors for this chart are obtained using the mean and standard deviation of the sample 
distribution of s which are fi s =c 2 <7 and <r s =C 3 a with: 


c 3 =<r 


n— 1 
y n 


-c 


2 • 


We have, then, Z?2 =c 2 + 3c3 and fii=C 2 - 3c3. 


(409) 

(410) 
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Given: s Central Line: s 


UCL: B 4 s 
LCL: Bj s 

where B 4 =B 2 /c 2 =(c 2 + 3c3)/c2 and Bj=B\/c 2 =(c 2 - 3cf)lc2- (411) 


J. Reliability and Life Testing 

In life testing , the random variable under consideration is the time-to-failure of a component. In 
fatigue studies, the time-to-failure is measured in the number of cycles to failure and is, therefore, a 
discrete random variable. 

The probability density /(r) associated with the time T to failure is called the failure-time density 
or life distribution. The probability that the component will fail in the time interval (0, t ) is given by the 
cumulative distribution: 


t 

P(T<t) = F(t) = \f(x)dx (412) 

0 

which, in this context, is called the unreliability function. The complementary cumulative distribution 
defines the probability that the component functions longer than time t or survives at least to time t. It is 
given by: 


P(T>t)=R{t)= l-F(f) . (413) 

This is called the reliability function by engineers and survivorship function by actuaries. 

NOTE: From a mathematical viewpoint, reliability is by definition a probability, therefore a number 
between 0 and 1 . 

By its very nature reliability theory has its focus in aerospace engineering on the tail area of a 
distribution, which defines low risk or high reliability systems. Thus the difference among different 
reliability models becomes significant only in the tails of the distribution where actual observations are 
sparse because of limited experimental data. In order to be able to discriminate between competing life 
distribution models, reliability engineers resort to a concept which attempts to differentiate among 
distributions because of physical failure mechanisms, experience, and/or intuition. Such a concept is the 
failure rate function which originated in actuary theory and mortality statistics where it is called the force 
of mortality. It is also indispensable in the analysis of extreme values where it is known as the intensity 
function. 

To define this function, we first determine the conditional probability that a component of “age” t 
will fail in the subsequent interval (r, t + At) by defining the following events: 

A=lifetime t <T<t + At (414) 

fi=lifetime T > t. (415) 
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Notice that the event B is the condition that no failure has occurred up to time f; i.e., the component has 
age t. In terms of these events the conditional probability is, therefore: 


P[A I B]= -P[t<T <t+At I T>t]= P[(t - T - p t f + T A ^ (T ~ t) U AG(t) 


(416) 


The numerator in the above equation is the probability that T lies between t and t+At and the denominator 
is the reliability function R(t). The equation can, therefore, be rewritten as: 


AG(t)= 


F(t+At)-F(t) 

R(t) 


(417) 


The failure rate is defined as the probability of failure in the interval per unit time given that no 
failure has occurred before time t. Accordingly, we divide the above equation by At and then take the limit 
At — » 0 to arrive at the failure rate function : 


Z(t)=F'( t)/R( t) = f(t)/R( t) = f(t)/{l- F(t )} . (418) 

This function is also called hazard function, hazard rate, or simply hazard. It is, in a sense, the 
propensity of failure as a function of age. The failure rate function contains the same information as the 
failure time distribution but it is better suited to formulate reliability models for given sets of experimental 
data. 


The difference between the failure time density fit) and the failure rate function z(t) can be 
illustrated by comparing it to human mortality. There /(f) dt is the probability that a newborn person will 

die between age t and t + At, whereas z(t) dt is the probability that a person of age t will die in the 
subsequent time interval At. 

Some statisticians have been investigating the reciprocal l/z(t) of the failure rate function, which is 
called Mill’s Ratio. It has no application in reliability studies. 

Two important relationships exist between the failure rate function and the reliability function. 
They can be obtained by the following steps. First using F(t)= 1 -R(t) we obtain by differentiation that 
F'(t) = -R'(t). 

Therefore: 


d[£nR(t)J 
Z[l> ~ R(t)~ dt 


(419) 


Solving this differential equation for R(t) and subsequently using the relationship f(t) =z(t ) R(t) 

yields: 


-\z{x)dx -\z(,x)dx 

R{t)-e 0 and f(t) = z(t)e 0 . (420) 
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The property that R (°°)=0 requires that the area under the failure rate curve be infinite. This is in 
distinction to the failure time density for which the area under its curve is one. Therefore the failure rate 
function is not a probability density. Whereas /(t) is the time rate of change of the unconditional failure 
probability, z(t ) is the time rate of change of the conditional failure probability given that the component 
has survived up to time t. 

EXAMPLE: The function f(t)=c e~ at does not quality as a failure rate function because the area under its 
curve is finite. Besides this requirement, the failure rate function has also to be nonnegative . Therefore we 
have the two properties of the hazard function: 

Z(t) > 0 for all t ( 421 > 

and 

S z ( t)dt = °° (422) 


1. Life Testing 

In life testing, a random sample of n components is selected from a lot and tested under specified 
operational conditions. Usually the life test is terminated before all units have failed. There are essentially 
two types of life tests. The first type is a test that is terminated after the first r failures have occurred (r < n) 
and the test time is random. Data obtained from such a test are called failure censored or type I censored 
data. The other more common type test is terminated after a predetermined test time. In this case, the 
number of failures is random. Such data are called time censored or type II censored data. Failure 
censoring is more commonly dealt with in the statistical literature because it is easier to treat 
mathematically. If the failed units are replaced by new ones, the life test is called a replacement test; 
otherwise, it is called a nonreplacement test. The unfailed units are variously called survivors, run-outs, 
removals, suspensions, or censored units. 

In recent years the Weibull distribution has become one of the most popular lifetime distributions. 
Because it can assume a wide variety of shapes, it is quite versatile. Although there exists a quick and 
highly visual graphical estimation technique, the maximum likelihood method is more accurate and lends 
itself easily to accommodate censored data. 

In the following, we present the maximum likelihood method for estimating the two parameters of 
the Weibull distribution and a BASIC computer program for solving the maximum likelihood equations. 

Let the life data consist of n units of which R units have failed at times r,- and S units have survived 
up to time ts . Assuming a Weibull distribution, the likelihood function is given by: 

L= Y[afitf~ X e~ m ^ Y[e~ at k - (423) 

i=l k = 1 

Taking the logarithm of the likelihood function called the log-likelihood function yields: 

R R S 

in L = R in a + R in ft + (/? - 1) ^Int^-a ^tf - a . (424) 

i'=l i=l Jt=l 
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The locations of the extreme values of the logarithm of a function are identical with those of the 
function itself. This can be shown by using the chain rule as follows: 


d£nF(x) = dinF(x) dFjx[ _ \ dF(x) 
dx dF dx F(x) dx 


(425) 


Therefore, if the function F(x) remains finite in the interval in which the local extreme values are located, 
the derivative of the logarithm of the function is zero at the same values of x for which the derivative of the 
function F(x) itself is zero. 


The derivatives with respect to the Weibull parameters a and /J are: 


and 


d In L 

da 


R 

a 


'z*f+ iff )=o 

\i=l K = 1 ) 


d in L R * 


= — + 


V P 


Y in tj - a\ 

i'=l 


X tf in tj + 
V(=l 


S 

I.K 

k= 


4 int^ = 0 . 


(426) 


(427) 


We can eliminate a from the first equation and insert it in the second equation to obtain a single 
nonlinear equation in /3. Thus, we have: 


a- 


R 


R „ S 
/=1 k = 1 


(428) 


and 


R 

P 


R 

+ 'L?nt •- 

i=l 


f R S 

/?l x *1 tn t { + X int k 


i=l 


k = 1 

~R~Z S~ 

i=i k=\ 


=o. 


(429) 


The second equation can be solved for (5 by an iterative method and then be used to solve for a. 

It is noteworthy that if there are no suspensions, it takes at least two failure times to be able to 
estimate the two Weibull parameters. However, if there are suspensions, only one failure time is sufficient 
to do this. On the other hand, one cannot obtain the parameters if there are only suspensions and no 
failures. That is why one often hears the statement that suspensions are a weak source of information. 

The following is a BASIC program using Newton’s method to calculate the maximum likelihood 
estimators for the Weibull parameters in an iterative manner. This method requires an initial guess for the 

shape parameter /? of the Weibull distribution. However, the initial values do not have to be too close to 
the actual parameter for the algorithm to converge. 
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WEIMLE 
DEFDBL A-Z 
INPUT"R=",R 
DIM STATIC R(R) 

FOR 1=1 TO R INPUT TR= , R(I):NEXT I 

130 :INPUT "S=’\S 
IF S=0 GOTO 140 
INPUT"TS=",TS 
'DIM STATIC S(S) 

TOR 1=1 TO S:INPUT TS= , S(I):NEXT I 


140 : INPUT "B=",B:J=0 

150 : H=B*.0001:GOSUB 500 

D=B:Y=F 

B=B+H:GOSUB 500 
B=D-H*Y/(F-Y):J=J+1 
IF ABS((D-B)/D)>. 000001 GOTO 150 
BEEP:BEEP 

PRINT "B=";B :PRINT ' J=" ;J 
A=R/(R 1 +S 1 ) : E= A A (- 1/B) 

PRINT" A=";A:PRINT"E="E 

END 


500 : R1=0:R2=0:R3=0 

FOR 1=1 TO R 

U=R(I) A B : V=LOG(R(I)) : W=U* V 
R1=R1+U:R2=R2+V:R3=R3+W :NEXT I 
S1=0:S3=0 
IF S=0 GOTO 580 

S 1=S*TS A B:S3=S 1 *LOG(TS):GOTO 580 
FOR 1=1 TO S 

U=S(I) A B:V=LOG(S(I)0:W=U*V 

S1=S1+U:S3=S3+W:NEXT I 

580 : F=R/B+R2-R*(R3+S3)/(R1+S1) 

RETURN 


2. No-Failure Weibull Model 

Sometimes when high reliability items are subjected to a time-censored test, no failures occur at 
all. This presents a real dilemma, because the traditional methods for estimating the parameters of a 
distribution cannot be applied. However, when it is assumed that the exponential model holds, one can 

determine an approximate lower (1 - a) confidence limit for the mean life of the component, which is 
given as: 


9T 

p> A 4 r with v= 2 degrees of freedom. 

xl 


(430) 
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Here, T is the fixed accumulated lifetime of the components tested and H>^ir cuts off a right- 

xl 

hand tail of the distribution with 2 degrees of freedom. This distribution is identical to the exponential 
distribution f{t)=2 e~ t/2 . Its critical value can be easily determined to be a. If there are n 

components, each having a lifetime /, on test, then the accumulated lifetime is: 


n 

T=lt r (431) 

i=l 

The lower (1 - a) confidence interval is then given by: 

n 

■¥ 

< 432 ) 

The upper confidence interval is, of course, infinity. 

PROBLEM: A sample of 300 units was placed on life test for 2,000 hours without a failure. Determine an 
approximately lower 95-percent confidence limit for its mean life. 

SOLUTION: 




(300)(2,000) _ 600,000 

-tn( 1-0.95) 2.995732274 


=200, 284 hours. 


(433) 


It is not known whether the approximation of the lower bound is conservative in the sense that the 
actual confidence is higher than the calculated one. In fact, some practicing statisticians advise against 
using this approximation at all. Admittedly, it is an act of desperation, but it appears to give reasonable 
results. 


The above method can be extended to the Weibull distribution by making use of the following 
relationship: If the failure time Thas a Weibull distribution with shape parameter /J and characteristic life 

T], then the random variable Y=T& has an exponential distribution with mean fi= l]P. To show this, we 
apply the Jacobian method, according to which we have: 


g(y)=f(t) 


dt 


dy 


= a[itP 


-1 


ptP~ 


Expressing the above equation in terms of the new variable y, we obtain: 


g(y)=ae °^with n y =±=r]P . 

Thus, if no failures have occurred in a time-censored test and the n components have lifetimes 
t\, t 2 ■ ■ . t n , the approximate lower confidence for ji y is: 


(434) 


(435) 
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(436) 


M 


— riP >J=1 


H v =n H 2 


-£/ia ' 


From this we obtain the lower confidence limit for the characteristic life: 

1 


f 


P 


ri> 


v 


-Ina . 


\P 


(437) 


Notice that the shape parameter (5 is still unknown and cannot be determined from the existing no-failure 
data set. To overcome this obstacle, it has been advocated to use historical failure data from similar 
components or from engineering knowledge of the failure physics under consideration to supply the 
missing Weibull shape parameter. This method is known as Weibayes analysis because it is regarded as 
an informal Bayesian procedure. 

Sometimes the denominator of the above equation is set equal to one. This corresponds to a 
(1 -e~ l ) or 63.2-percent confidence interval for the lower bound. It has been said that for this condition 
“the first failure is imminent.” This statement must surely rank among the most bizarre ever made in life 
testing analysis. 


K. Error Propagation Law 

Instead of aspiring to calculate the probability distribution of an engineering system response or 
performance from the distribution of the component variables (life drivers), it is frequently adequate to 
generate some of its lower order moments. It is known, for instance, that the knowledge of the first four 
moments of a distribution enables one to establish a Pearson-type distribution or a generalized lambda 
distribution from which percentiles of the distribution can be estimated. Since these two families of 
distributions can accommodate a wide variety of distributions, the technique is often quite effective and the 
result is usually better than expected. Generating approximate system performance moments by this 
method is generally known as statistical error propagation. Sometimes it is referred to as the Delta 
Method or the Root-Sum-Square (RSS) Method. 

In order to apply this technique, the functional relationship between the system output 
(performance) and system parameters (life drivers) must be explicitly known or at least be approximated. 
The essence of the technique consists in expanding the functional relation in a multivariate Taylor series 
about the design point or expected value of the system parameters and retaining only lower order terms. 
The goodness of the approximation depends strongly on the magnitude of the coefficient of variation of the 
system parameter distributions; the smaller the coefficient of variation, the better is the approximation. 

The analysis is actually rather straightforward, albeit tedious, for moments higher than the second. 
A software package called Second Order Error Propagation (SOERP) relieves one of this drudgery. 
Despite its misleading acronym, it performs higher order error propagation analysis for complicated 
engineering systems. The following exposition is limited to calculating only the first two moments, namely 
the mean and variance, of a system output. 
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Let the performance of a system in terms of its system parameters jq, X 2 . . . x n be given by: 

z=f(x i, *2 •••*«) • (438) 



Let H=(fl i, H2 - ■ Hn) and 0=(0\> 02 • • • 0n) denote the vector of the mean and standard deviation 
of the system parameters, respectively. 

The Taylor series expansion of the function z about the point [i 2 ...//„) can be 

symbolically defined by: 


Z = f(x v x 2 ,...x n )= I -7 
n = 0 "• 


Ax, — hAx-, 


in 


>5^ 2 C^ 


. .Ax, 


n dx„ 


n M’V 2 -V n 


f(x v x 2 ...x n ) (439) 


where Ax,- (x, ■-/!,•) and the partial derivatives are evaluated at the design point fj. 2 , 

Retaining only terms up to second order we have: 


z f(H v /i 2 ...fi n )+ ( x i (*i Vi)(Xj flj). 

ft f. I 


d 2 f 


(440) 


1. The Mean of the System Performance 

To calculate the mean of the system performance we take the expectation of the above equation and 
observing that E((Xj-fii))=0, we obtain: 


E(z)=H z =f(H)+±i 


i=l j=\\d*i dx j ) 


E[(x i -fj, i )(x j -Hj)] . 


(441) 


Deleting the second-order term, we get the first-order approximation for the mean system performance: 

E(z)=fl z =f(fi) . (442) 

2. The Variance of the System Performance 


The calculation of the variance is more complex than the one of the mean. In fact, the calculations 
become substantially more lengthy with increasing order of the moments. To begin with, we calculate the 
variance for the case of two system parameters x and y. This is not a real loss in generality. Once the result 
for this case has been obtained, the extension to more than two variables can be rather easily recognized. 
Let the system performance function be given by: 


z=f(x,y ) . (443) 

Furthermore, let us denote the variance of a random variable X by V(X)=E(X — E(X))^=E(X - 
Hxp-=cft and the partial derivative of a function with respect to a variable, say x, by f x . 

(a) First-Order Taylor Series Approximation. Here we take only the first terms of the above 
Taylor series. The first-order approximation of the variance of the system performance is then obtained by: 
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(444) 


V(z)=E(z-H z ) 2 =E[fx(x-n x )+fy(y-Li y )f 

= fx °x + fy a l +2 fxfy °xy 

where o xy is the covariance of x and y. The partial derivatives f x and f y are sometimes called sensitivity 
coefficients. 


Recalling that the covariance is related to the correlation coefficient by <7xy~P &x Oy, we obtain the 
standard deviation of the system performance by taking the square root of the above equation: 


C Z=pl °l + fy a y +2 fxfyP°x a y • ( 445 ) 

Extension of this expression to more than two system parameters is straightforward. An important, 

because often encountered, situation arises if the system parameters are independent (p= 0). For this 
special case the standard deviation of the system performance simplifies to: 


(j = f 2 cj-2 _|_ rz fji 

z \ J x x J y y 


r 2 „2 


(446) 


This formula is the reason why this method is also called Root-Sum-Square Method. If correlation 
is suspected but unknown, one can use the upper limit of the standard deviation which is given by: 


O z <\f x \a x +\f y \a y . (447) 

This corresponds to a “worst-on-worsf ’ type situation. 

Before proceeding to the second-order approximation, we examine some special cases of practical 
importance. 

Multi-Output Linear System y=A x. Let the system performance function be a vector function of 
the following form: 


y=f(x) ■ (448) 

where the system performance (output) vector y is an (m x 1) vector and the system parameter vector x be 

an (n x 1) vector. We expand now the system performance function in a Taylor series. Without loss of 
generality we set the value of the performance function at the design point equal to zero and obtain the 
linearized system performance function as: 


y=A x (449) 

where A is an ( m x n) matrix. It is known as the Jacobian matrix or system sensitivity matrix and its 
elements are the first partial derivatives of the system performance vector function y=f(x). 

To calculate the first-order approximation of the system performance variance, we introduce the 
covariance matrix which is defined as: 
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(450) 


G\ (7\2---<7\n 
On c\...Cl n 

<7 In (72 n ••• On . 

and is symmetric. 

The covariance matrix of the system performance is then obtained by taking the following steps: 
Cov (y) = E [(x-tty) ix-Uy^E [XX 7 ) ~ E[x\ E\x T 1 

= E [A xx t A t ] - E [Ax] E [x T A?] . (451) 

Continuing from the previous page: 

Cov (y) = AE[xx T ]A T -AE[x]E[x T ]A T 

= A [E (x x T ) - E (y) E (y 7 )] A T (452) 

and finally 

Cov (y) = A [Cov(y)] AT . (453) 


As an example: 

yi = a\ X] + ci 2 X 2 (454) 

X2 = b\x\ + t >2 X 2 . (455) 

We assume independent system parameters; i.e., p= 0. 


or 


Cov (y) 


Cov (y) 


a l “2 

1 

Q 

o 

1 

~ a l 

V 

b 2_ 

L° a i\ 

“2 

hi 

a \ °2 

a \(7\ fcjCTj 2 


b 1 b 2_ 

_ a 2°2 b 2°2 . 


( 

Q 

+ a l<- 7 i) 



(456) 


(457) 


(458) 


\[ a \ b \ a \ +a 2 b 2 a 2) { b \°\+ b 2 a 2) J 

NOTE: Although the system parameters were assumed to be independent, the system performance output 
variables become dependent because of the “cross-coupling” (off-diagonal) terms of the sensitivity matrix. 
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Power Function z=x m y n . Let us denote the coefficient of variation as ri=a/fi. The partial 
derivatives of the system parameters evaluated at the design point |l=(|l x , p.y) are: 

fx=m (^r- 1 <Jl y ) n and f y =n QiJ" (/ty)"- 1 . (459) 

Assuming again independence of the system parameters, the system performance variance is then 
given by: 

0*=[m(n x ) m - l (H y ) n ] 0 2 +\n(n x ) m (ll y ) n ~ 1 ] crj. (460) 

Dividing both sides by n 2 =^(/J. x ) m (li y ) n j we obtain an expression for the coefficient of 
variation for the system performance which is: 

ri z = ^m 2 Tll+n 2 T]2 . (461) 

We see that the coefficient of variation of the system performance can be expressed in terms of the 
coefficient of variation of the system parameters and that it is magnified by the power of the system 
parameters. 

For example, the volume of a sphere is V=4/3 tt/? 3 where R is its radius. Therefore a 10-percent error in 
the radius will show up as a 30-percent error in the volume of the sphere. 

(b) Second-Order Taylor Series Approximation. Let us first consider the case of only one 
system parameter. The Taylor series expansion about the mean j u x is: 

* = /og + /, <*-/g ■ + 'A <■* ■ < 462 > 

The variance is best calculated by using the relationships 

V(a + jc)=V(x) a=constant (463) 

and 

V{x)=E (x 2 )-{E {x)) 2 . (464) 

The mean is simply given by the expectation of the Taylor series, which is: 

V z =f(V x )+±fxx a l- ( 465 > 

Therefore, the variance of the system performance is 

(466) 

= f 2 E(x-^l x ) 1 +\fl c Etx-^f+fJn E(x-ll x ) 2 > -\fl c E 2 (x-h x ) 2 
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( 467 ) 


°Z =f?° 2 x + \fxx^ + fxfxx»3-\fl a x 
where jU 3 =skewness and ^ 4 =kurtosis of the system parameter. 

For example: z= x 2 - 

We assume the system parameter to be normally distributed with /J. x =0 and cr^=l . We have the 

derivatives f x =2x,f xx =2 and the moments /J x =0, <r^=l , Hi=0, and fi4=3. From this we obtain the mean 
and variance of the system performance as 


Hz=l and o 2 x = 2 . (468) 

It can be shown that these moments agree with the exact moments of the variable z, which in this case has 
a x 2 distribution with 1 degree of freedom. 

The case of several parameters becomes more complicated and will not be further pursued. For 
uncorrelated system parameters, the resulting expressions for the system performance moments, retaining 
only third-order terms in the final result, yield: 


Variance: 


= X 


n f \ 


* 

dX; 





dxf 

v * y 


H (*i) 


(469) 


n ( 


Skewness: 


* (x i> 


nfdf^ 


Kurtosis: /2 4 (z)=Z, 

i=ly (XX i 


/M*») + 6 X I 


' J 

i>j 


f^l 

2 

f 





a 1 (Xi)o 2 (xj) . 


(470) 


(471) 


All derivatives are again evaluated at the design point of the system. 
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